To run: Beforehand:
module load pandoc
In R:
setwd("~/shared-gandalm/brain_CTP/Scripts/methylation/analysis/combine_methods/ctp_validation")
rmarkdown::render("ctp_validation_ROSMAP.Rmd", "html_document")
library(tidyverse)
library(rmarkdown) # You need this library to run this template.
library(epuRate) # Install with remotes::install_github("holtzy/epuRate", force=TRUE) - use this instead of devtools::install_github()
library(data.table)
library(DT)
library(tidyverse)
library(dplyr)
library(reshape2)
library(uwot) # for umap
library(methylCC)
library(lme4)
library(compositions)
library(zCompositions)
library(mixOmics)
source("~/project-gandalm/software/ANCOM-v2.1/scripts/ancom_v2.1.R")
library(ALDEx2)
library(factoextra)
library(ggplot2)
library(RColorBrewer)
library(viridis)
library(ggpubr)
library(GGally)
library(lattice)
library(gridExtra)
# BiocManager::install("M3C")
# library(M3C)
#filen <- "ROSMAP_arrayMethylation_imputed"
#filen <- "ROSMAP_ComBatbatchplate_neg0"
#filen <- "ROSMAP_ComBatplate_neg0"
#filen <- "ROSMAP_SNMplate_neg0"
filen <- "ROSMAP"
# Cell-type number
cellnum <- 9
cellnum <- 7
write_out = TRUE
out_dir <- "~/shared-gandalm/brain_CTP/Data/methylation/ROSMAP/analysis"
# Read in reference object
# - use chisq statistic, 5 cell-types, bonferroni p<0.001, marker contributing >=50% of chisq stat, take top 200
# ref_dir <- "~/shared-gandalm/brain_CTP/Data/reference_cell_profile/Luo2020/Luo2020_chisq_dmr_ilmn450kepic_aggto100bp_cell7_bonf001_contr50pc_top200.rds"
# ref_dir <- "~/shared-gandalm/brain_CTP/Data/reference_cell_profile/Luo2020/Luo2020_chisq_dmr_cell5_bonf001_contr50pc_top200.rds"
#ref_dir <- "~/shared-gandalm/brain_CTP/Data/reference_cell_profile/Luo2020/Luo2020_chisq_dmr_ilmn450kepic_aggto100bp_cell7_bonf001_contr50pc_top200.rds"
#ref <- readRDS(ref_dir)
# Pheno
pheno_dir <- "~/shared-gandalm/brain_CTP/Data/pheno/ROSMAP/pheno_ROSMAP.txt"
# methylCC
# - use Luo et al. reference with DMRs identified from 450K overlap
#saveRDS(methylcc_out, paste(out_dir, "/", filen, "_aut_mask_methylcc_out_n200_p1e-6.rds", sep = ""))
# - use 450K dlpfc guintivano reference (+/- randomisation of sample order)
#saveRDS(methylcc_out, paste(out_dir, "/", filen, "_aut_mask_methylcc_out_n200_dlpfc_450k_guintivano_hpc.rds", sep = ""))
# - use Luo et al. reference with DMRs pre-identified by Luo et al.
#saveRDS(methylcc_out, paste(out_dir, "/", filen, "_aut_mask_predmr_ilmn450k_celltype_na5_methylcc_dmr_n200_p1e_1.rds", sep = ""))
# - use extreme_methylation_sites.R output
#saveRDS(methylcc_out, paste(out_dir, "/", filen, "_aut_mask_extreme17_dmr_low0.3_high0.7.rds", sep = ""))
# - use chisq statistic, 5 cell-types, bonferroni p<0.001, marker contributing >=50% of chisq stat, take top 200
#saveRDS(methylcc_out, paste(out_dir, "/", filen, "_aut_mask_chisq_dmr_cell5_bonf001_contr50pc_top200.rds", sep = ""))
#methylcc_dir <- paste(out_dir, "/", filen, "_aut_mask_chisq_dmr_cell5_bonf001_contr50pc_top200.rds", sep = "")
filen0 <- "ROSMAP"
#methylcc_dir <- paste(out_dir, "/", filen0, "_aut_mask_chisq_dmr_cell5_bonf001_contr50pc_top200_aggpostmcc.txt", sep = "")
#methylcc_dir <- paste(out_dir, "/", filen0, "_aut_mask_chisq_dmr_ilmn450kepic_aggto100bp_c10_cell7_bonf001_contr50pc_top200_aggpostmcc.txt", sep = "")
if (cellnum == 9) {
methylcc_dir <- paste(out_dir, "/", filen0, "_aut_mask_extremes_dmr_ilmn450kepic_aggto100bp_c10_cell9_split6040all_aggpostmcc.txt", sep = "")
} else if (cellnum == 7) {
methylcc_dir <- paste(out_dir, "/", filen0, "_aut_mask_extremes_dmr_ilmn450kepic_aggto100bp_c10_cell7_split6040all_aggpostmcc.txt", sep = "")
}
# methylCC with NeuN+/- data (Guintivano DLPFC)
methylcc_neun_dir <- paste("~/shared-gandalm/brain_CTP/Data/methylation/ROSMAP/analysis/", filen, "_aut_mask_methylcc_out_n200_dlpfc_450k_guintivano_hpc.rds", sep = "")
# eBayes + Houseman
houseman_dir <- paste("~/shared-gandalm/brain_CTP/Data/methylation/ROSMAP/analysis/", filen, "_aut_mask_t_ref450K_toast_houseman_ls.rds", sep = "")
# sequencing with extremes + Houseman
houseman_seq_dir <- "~/shared-gandalm/brain_CTP/Data/methylation/ROSMAP/analysis/ROSMAP_aut_mask_t_Luo2020_extremes_dmr_ilmn450kepic_aggto100bp_c10_cell7_split6040all_houseman_ls.rds"
# celfie
celfie_dir <- "~/shared-gandalm/brain_CTP/Data/methylation/ROSMAP/analysis/celfie/ROSMAP_aut_mask_t_celfie.out/1_tissue_proportions.txt"
celfie_label_dir <- "~/shared-gandalm/brain_CTP/Data/methylation/ROSMAP/analysis/celfie/ROSMAP_aut_mask_t_celfie.in"
# ReFACToR
#refactor_dir <- "~/shared-gandalm/brain_CTP/Data/methylation/ASD_methylation_brain/analysis/KCL_R01MH094714_ASD_Illumina450K_PFC_k8_cov.refactor.components.txt"
# smartSV
filen <- "ROSMAP_ComBatplate_neg0"
smartsva_dir <- paste("~/shared-gandalm/brain_CTP/Data/methylation/ROSMAP/analysis/", filen, "_aut_mask_t_SmartSVA.rds", sep = "")
filen <- "ROSMAP"
# VAE/MethylNet
vae_dir <- "~/shared-gandalm/brain_CTP/Data/integration/methylnet/combined/full_output_latent.csv"
# Pheno
pheno <- read.delim(pheno_dir, header = T, as.is = T)
colnames(pheno)[which(colnames(pheno) == "individualID")] <- "IID"
colnames(pheno)[which(colnames(pheno) == "age_death")] <- "age"
pheno$age[which(pheno$age == "90+")] <- 90
pheno$age <- as.numeric(pheno$age)
colnames(pheno)[which(colnames(pheno) == "msex")] <- "sex"
pheno$sex <- ifelse(pheno$sex == 1, "M", "F")
pheno$batch <- as.factor(pheno$batch)
pheno$group <- ifelse(pheno$cogdx %in% 4:5, "AZD", "CTL")
# methylCC
methylcc_ctp <- read.delim(methylcc_dir)
ind <- grep("NonN", colnames(methylcc_ctp))
foo <- colnames(methylcc_ctp)[ind]
foo2 <- unlist(lapply(strsplit(gsub("NonN_", "", foo), "_"), function(x) x[1]))
colnames(methylcc_ctp)[ind] <- foo2
colnames(methylcc_ctp)[2:ncol(methylcc_ctp)] <- paste("mcc_", colnames(methylcc_ctp)[2:ncol(methylcc_ctp)], sep = "")
# methylcc <- readRDS(methylcc_dir)
# methylcc_ctp <- cell_counts(methylcc)
# ctp_sum <- methylcc@summary
# ctp_input <- ctp_sum$input
# colSums(ctp_input$zmat)
# Would use this to mimic Supp Fig 1
# colSums(ref$zmat) # this is the initial starting proportions
# methylCC on Guintivano DLPFC reference (NeuN+/-)
# methylcc_neun <- readRDS(methylcc_neun_dir)
# methylcc_neun_ctp <- cell_counts(methylcc_neun)
# colnames(methylcc_neun_ctp) <- c("mcc_NeuN_neg", "mcc_NeuN_pos")
# ctp_neun_sum <- methylcc_neun@summary
# ctp_neun_input <- ctp_neun_sum$input
# colSums(ctp_neun_input$zmat)
# eBayes + Houseman
houseman.ls <- readRDS(houseman_dir)
# - coefs_sub
coefs_sub <- houseman.ls$coefs_sub
colnames(coefs_sub)[2:ncol(coefs_sub)] <- paste("h_", colnames(coefs_sub)[2:ncol(coefs_sub)], sep = "")
colnames(coefs_sub)[ncol(coefs_sub)] <- "h_error_sub"
coefs_brain <- houseman.ls$coefs_brain
colnames(coefs_brain)[2:ncol(coefs_brain)] <- paste("h_", colnames(coefs_brain)[2:ncol(coefs_brain)], sep = "")
colnames(coefs_brain)[ncol(coefs_brain)] <- "h_error_brain"
# sequencing + Houseman
houseman_seq.ls <- readRDS(houseman_seq_dir)
houseman_seq <- houseman_seq.ls[[1]]
ind <- grep("NonN", colnames(houseman_seq))
foo <- colnames(houseman_seq)[ind]
foo2 <- unlist(lapply(strsplit(gsub("NonN_", "", foo), "_"), function(x) x[1]))
colnames(houseman_seq)[ind] <- foo2
colnames(houseman_seq)[2:ncol(houseman_seq)] <- paste("hseq_", colnames(houseman_seq)[2:ncol(houseman_seq)], sep = "")
# CelFIE
celfie <- fread(celfie_dir)
colnames(celfie) <- c("IID", "celfie_Exc", "celfie_Inh", "celfie_Astro", "celfie_Endo", "celfie_Micro", "celfie_Oligo", "celfie_OPC", "unknown1")
#celfie_label <- fread(celfie_label_dir)
# ReFACToR
# refactor <- read.table(refactor_dir, header = T, as.is = T)
# colnames(refactor)[1] <- "IID"
# smartSVA
smartsva <- readRDS(smartsva_dir)
smartsva <- data.frame(IID = rownames(smartsva), smartsva)
if (ncol(smartsva) > 11) {
smartsva <- smartsva[,c(1:11)]
}
# VAE, MethylNet
vae <- read.csv(vae_dir)
colnames(vae) <- c("IID", paste("VAEe", 0:9, sep = ""))
#ctp_pheno <- inner_join(data.frame(IID = rownames(methylcc_ctp), methylcc_ctp), pheno, by = "IID")
methylcc_ctp.tmp <- methylcc_ctp
methylcc_ctp.tmp$mcc_Glial <- rowSums(methylcc_ctp.tmp[,-c(grep("IID", colnames(methylcc_ctp.tmp)), grep("Exc", colnames(methylcc_ctp.tmp)), grep("Inh", colnames(methylcc_ctp.tmp)))])
methylcc_ctp.tmp$mcc_Neuron <- rowSums(methylcc_ctp.tmp[,c(grep("Exc", colnames(methylcc_ctp.tmp)), grep("Inh", colnames(methylcc_ctp.tmp)))])
ctp_pheno <- inner_join(methylcc_ctp.tmp, pheno, by = "IID")
#ctp_pheno <- inner_join(ctp_pheno, data.frame(IID = rownames(methylcc_neun_ctp), methylcc_neun_ctp), by = "IID")
ctp_pheno <- inner_join(ctp_pheno, smartsva, by = "IID")
#ctp_pheno <- inner_join(ctp_pheno, houseman0, by = "IID")
# Houseman with array reference
ctp_pheno <- inner_join(ctp_pheno, coefs_sub, by = "IID")
ctp_pheno <- inner_join(ctp_pheno, coefs_brain, by = "IID")
ctp_pheno$h_Glial <- rowSums(ctp_pheno[,c("h_ASTRO", "h_OPC", "h_OLIGO", "h_ENDO", "h_MONO")])
ctp_pheno$h_Neuron <- rowSums(ctp_pheno[,c("h_GABA", "h_GLU")])
# Houseman with sequencing reference
ctp_pheno <- inner_join(ctp_pheno, houseman_seq, by = "IID")
ctp_pheno$hseq_Glial <- rowSums(ctp_pheno[,c("hseq_Astro", "hseq_Endo", "hseq_Micro", "hseq_Oligo", "hseq_OPC")])
ctp_pheno$hseq_Neuron <- rowSums(ctp_pheno[,c("hseq_Exc", "hseq_Inh")])
# CelFIE
ctp_pheno <- inner_join(ctp_pheno, celfie, by = c("IID"))
ctp_pheno$celfie_Glial <- rowSums(ctp_pheno[,c("celfie_Astro", "celfie_Endo", "celfie_Micro", "celfie_Oligo", "celfie_OPC")])
ctp_pheno$celfie_Neuron <- rowSums(ctp_pheno[,c("celfie_Exc", "celfie_Inh")])
ctp_pheno <- inner_join(ctp_pheno, vae, by = "IID")
if (write_out == TRUE) {
write.table(ctp_pheno, paste(out_dir, "/", filen, "_ctp_pheno.txt", sep = ""), col.names = T, row.names = F, quote = F, sep = "\t")
}
keep_cols <- c("IID", "group", "age", "sex", "batch", "Sample_Plate", "race", "pmi")
ctp_cols <- c(colnames(ctp_pheno)[grep("hseq", colnames(ctp_pheno))])
ctp_cols <- ctp_cols[-which(ctp_cols == "hseq_error")]
level2_cols <- ctp_cols[-which(ctp_cols %in% c("hseq_Glial", "hseq_Neuron"))]
glial_rmoligo_n <- colnames(ctp_pheno)[grep("hseq", colnames(ctp_pheno))]
glial_rmoligo_n <- glial_rmoligo_n[-which(glial_rmoligo_n %in% c("hseq_Exc", "hseq_Inh", "hseq_Oligo", "hseq_error", "hseq_Glial", "hseq_Neuron"))]
neuron_n <- c("hseq_Exc", "hseq_Inh")
# Batch
batch <- ctp_pheno %>% group_by(group, batch) %>% dplyr::summarise(n=n()) %>% spread(group, n)
## `summarise()` has grouped output by 'group'. You can override using the `.groups` argument.
batch
## # A tibble: 2 × 3
## batch AZD CTL
## <fct> <int> <int>
## 1 0 100 164
## 2 1 199 252
chisq.test(batch %>% dplyr::select("AZD", "CTL"))
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: batch %>% dplyr::select("AZD", "CTL")
## X-squared = 2.419, df = 1, p-value = 0.1199
# PMI
summary(lm(pmi ~ group, data = ctp_pheno))
##
## Call:
## lm(formula = pmi ~ group, data = ctp_pheno)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.802 -3.052 -1.752 0.931 77.281
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.9545 0.3327 20.902 <2e-16 ***
## groupCTL 0.8478 0.4362 1.944 0.0523 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.753 on 713 degrees of freedom
## Multiple R-squared: 0.005271, Adjusted R-squared: 0.003876
## F-statistic: 3.778 on 1 and 713 DF, p-value: 0.05233
ggplot(ctp_pheno, aes(x = pmi, colour = group, fill = group)) + geom_histogram() + facet_wrap(~group)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#chisq.test(batch %>% dplyr::select("ASD", "CTL"))
# Race
race <- ctp_pheno %>% group_by(group, race) %>% dplyr::summarise(n=n()) %>% spread(group, n)
## `summarise()` has grouped output by 'group'. You can override using the `.groups` argument.
race
## # A tibble: 4 × 3
## race AZD CTL
## <int> <int> <int>
## 1 1 292 406
## 2 2 7 8
## 3 3 NA 1
## 4 7 NA 1
#chisq.test(batch %>% dplyr::select("ASD", "CTL"))
# Sex
sex <- ctp_pheno %>% group_by(group, sex) %>% dplyr::summarise(n=n()) %>% spread(group, n)
## `summarise()` has grouped output by 'group'. You can override using the `.groups` argument.
sex
## # A tibble: 2 × 3
## sex AZD CTL
## <chr> <int> <int>
## 1 F 199 255
## 2 M 100 161
chisq.test(sex %>% filter(sex %in% c("M", "F")) %>% dplyr::select("AZD", "CTL"))
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: sex %>% filter(sex %in% c("M", "F")) %>% dplyr::select("AZD", "CTL")
## X-squared = 1.8537, df = 1, p-value = 0.1734
# Age
summary(lm(age ~ group, data = ctp_pheno))
##
## Call:
## lm(formula = age ~ group, data = ctp_pheno)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.225 -2.246 1.983 2.553 4.782
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 87.8331 0.2650 331.499 < 2e-16 ***
## groupCTL -2.6152 0.3474 -7.529 1.55e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.582 on 713 degrees of freedom
## Multiple R-squared: 0.07364, Adjusted R-squared: 0.07234
## F-statistic: 56.68 on 1 and 713 DF, p-value: 1.552e-13
ggplot(ctp_pheno, aes(x = age, colour = group, fill = group)) + geom_histogram() + facet_wrap(~group)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# Education
summary(lm(educ ~ group, data = ctp_pheno))
##
## Call:
## lm(formula = educ ~ group, data = ctp_pheno)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.3712 -2.4832 -0.3712 2.5168 11.6288
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 16.3712 0.2071 79.067 <2e-16 ***
## groupCTL 0.1119 0.2715 0.412 0.68
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.58 on 713 degrees of freedom
## Multiple R-squared: 0.0002384, Adjusted R-squared: -0.001164
## F-statistic: 0.17 on 1 and 713 DF, p-value: 0.6802
ggplot(ctp_pheno, aes(x = educ, colour = group, fill = group)) + geom_histogram() + facet_wrap(~group)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# Braak stage
braaksc <- ctp_pheno %>% group_by(group, braaksc) %>% dplyr::summarise(n=n()) %>% spread(group, n)
## `summarise()` has grouped output by 'group'. You can override using the `.groups` argument.
braaksc
## # A tibble: 7 × 3
## braaksc AZD CTL
## <int> <int> <int>
## 1 0 NA 9
## 2 1 12 50
## 3 2 11 61
## 4 3 65 147
## 5 4 82 119
## 6 5 122 30
## 7 6 7 NA
# CERAD score (neuritic plaques)
ceradsc <- ctp_pheno %>% group_by(group, ceradsc) %>% dplyr::summarise(n=n()) %>% spread(group, n)
## `summarise()` has grouped output by 'group'. You can override using the `.groups` argument.
ceradsc
## # A tibble: 4 × 3
## ceradsc AZD CTL
## <int> <int> <int>
## 1 1 149 68
## 2 2 103 136
## 3 3 20 55
## 4 4 27 157
tmp.long <- ctp_pheno %>% dplyr::select(c("IID", "batch", "pmi", "Sample_Plate", "hseq_Glial", "sSV1", "sSV2", "sSV3", "sSV4")) %>% melt(id.vars = c("IID", "batch", "pmi", "Sample_Plate", "hseq_Glial"), variable.name = "sSV_num", value.name = "sSV")
ggplot(tmp.long, aes(x = batch, y = sSV)) + geom_boxplot() + facet_wrap(~ sSV_num, scales = "free_y")
ggplot(tmp.long, aes(x = pmi, y = sSV, colour = batch)) + geom_point() + facet_wrap(~ sSV_num)
ggplot(tmp.long, aes(x = Sample_Plate, y = sSV, colour = batch)) + geom_boxplot() + facet_wrap(~ sSV_num)
ggplot(tmp.long, aes(x = hseq_Glial, y = sSV, colour = batch)) + geom_point() + geom_smooth(method = "lm") + facet_wrap(~ sSV_num)
## `geom_smooth()` using formula 'y ~ x'
tmp.long <- ctp_pheno %>% dplyr::select(c("IID", "batch", "pmi", "age", "sex", "group", "Sample_Plate", "hseq_Glial", "VAEe0", "VAEe1", "VAEe2", "VAEe3", "VAEe4", "VAEe5", "VAEe6", "VAEe7", "VAEe8", "VAEe9")) %>% melt(id.vars = c("IID", "batch", "pmi", "age", "sex", "group", "Sample_Plate", "hseq_Glial"), variable.name = "VAEe", value.name = "embedding")
ggplot(tmp.long, aes(x = batch, y = embedding)) + geom_boxplot() + facet_wrap(~ VAEe, scales = "free_y")
ggplot(tmp.long, aes(x = age, y = embedding, colour = batch)) + geom_point() + facet_wrap(~ VAEe)
ggplot(tmp.long, aes(x = pmi, y = embedding, colour = batch)) + geom_point() + facet_wrap(~ VAEe)
ggplot(tmp.long, aes(x = Sample_Plate, y = embedding, colour = batch)) + geom_boxplot() + facet_wrap(~ VAEe)
ggplot(tmp.long, aes(x = sex, y = embedding)) + geom_boxplot() + facet_wrap(~ VAEe)
ggplot(tmp.long, aes(x = group, y = embedding)) + geom_boxplot() + facet_wrap(~ VAEe, scales = "free_y")
ggplot(tmp.long, aes(x = hseq_Glial, y = embedding, colour = batch)) + geom_point() + geom_smooth(method = "lm") + facet_wrap(~ VAEe)
## `geom_smooth()` using formula 'y ~ x'
Add analysis info on
houseman_seq.long <- melt(houseman_seq, variable.name = "celltype", value.name = "seq_houseman_CTP")
## Using IID as id variables
houseman_seq.gg <- ggplot(houseman_seq.long, aes(x = celltype, y = seq_houseman_CTP, colour = celltype)) + geom_boxplot() + theme(legend.position = "none", axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) + ylim(c(0,1))
methylcc_ctp.long <- melt(methylcc_ctp, variable.name = "celltype", value.name = "methylcc_CTP")
## Using IID as id variables
methylcc_ctp.gg <- ggplot(methylcc_ctp.long, aes(x = celltype, y = methylcc_CTP, colour = celltype)) + geom_boxplot() + theme(legend.position = "none", axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) + ylim(c(0,1))
#+ theme(legend.position = "none", axis.text.x = element_text(angle = 45, hjust = 1))
coefs_sub.long <- melt(coefs_sub, variable.name = "celltype", value.name = "eBayes_houseman_CTP") %>% mutate(MatchCell = match(celltype, c("h_GLU", "h_GABA", "h_ASTRO", "h_ENDO", "h_MONO", "h_OLIGO", "h_OPC")))
## Using IID as id variables
coefs_sub.gg <- coefs_sub.long %>% mutate(celltype = fct_reorder(celltype, MatchCell)) %>% ggplot(aes(x = celltype, y = eBayes_houseman_CTP, colour = celltype)) + geom_boxplot() + theme(legend.position = "none", axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) + ylim(c(0,1))
coefs_brain.long <- melt(coefs_brain, variable.name = "celltype", value.name = "eBayes_houseman_CTP")
## Using IID as id variables
coefs_brain.gg <- ggplot(coefs_brain.long, aes(x = celltype, y = eBayes_houseman_CTP, colour = celltype)) + geom_boxplot() + theme(legend.position = "none", axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) + ylim(c(0,1))
smartsva.long <- melt(smartsva, variable.name = "celltype", value.name = "smartsva_CTP")
## Using IID as id variables
smartsva.gg <- ggplot(smartsva.long, aes(x = celltype, y = smartsva_CTP, colour = celltype)) + geom_boxplot() + theme(legend.position = "none", axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
vae.long <- melt(vae, id.vars = c("IID"), variable.name = "celltype", value.name = "VAE_embed")
vae.gg <- ggplot(vae.long, aes(x = celltype, y = VAE_embed, colour = celltype)) + geom_boxplot() + theme(legend.position = "none", axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) + ylim(c(0,1))
celfie.long <- melt(celfie, id.vars = c("IID"), variable.name = "celltype", value.name = "CelFIE_CTP")
celfie.gg <- ggplot(celfie.long, aes(x = celltype, y = CelFIE_CTP, colour = celltype)) + geom_boxplot() + theme(legend.position = "none", axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) + ylim(c(0,1))
# Arrange into 1 plot
#ggarrange(houseman_seq.gg, methylcc_ctp.gg, coefs_sub.gg, coefs_brain.gg, jaffe.gg, celfie.gg, smartsva.gg, nrow = 3, ncol = 2)
ggarrange(houseman_seq.gg, methylcc_ctp.gg, coefs_sub.gg, coefs_brain.gg, celfie.gg, smartsva.gg, nrow = 3, ncol = 2)
## Warning: Removed 1 rows containing non-finite values (stat_boxplot).
# - boxplot
hseq_ctp_pheno.tmp <- ctp_pheno %>% dplyr::select(c(all_of(keep_cols), all_of(ctp_cols)))
if (write_out == TRUE) {
write.table(hseq_ctp_pheno.tmp, paste(out_dir, "/hseq_ctp_pheno.raw.txt", sep = ""), col.names = T, row.names = F, quote = F, sep = "\t")
}
hseq_ctp_pheno.long <- melt(hseq_ctp_pheno.tmp, measure.vars = c(all_of(ctp_cols)), variable.name = "celltype", value.name = "CTP")
hseq_ctp_pheno.long$Level <- ifelse(hseq_ctp_pheno.long$celltype %in% c("hseq_Neuron", "hseq_Glial"), "Level1", "Level2")
hseq_ctp_pheno.long$Region <- "PFC"
hseq_ctp_pheno.long$celltype <- gsub("hseq_", "", hseq_ctp_pheno.long$celltype)
hseq_ctp_pheno.long$celltype <- unlist(lapply(strsplit(hseq_ctp_pheno.long$celltype, "_"), function(x) x[1]))
hseq_ctp_pheno.long$MatchCell <- match(hseq_ctp_pheno.long$celltype, rev(c("Exc", "Inh", "Astro", "Micro", "Endo", "Oligo", "OPC", "Neuron", "Glial")))
hseq_ctp_pheno.long %>% mutate(celltype = fct_reorder(celltype, MatchCell)) %>%
ggplot(aes(x=celltype, y = CTP, fill = group)) +
geom_boxplot(position = position_dodge(.8),alpha=.8,outlier.size = 0)+
geom_point(position = position_jitterdodge(.1),alpha=.3) +
facet_grid(Level~Region, scales='free', space = 'free_y') + coord_flip() + theme_bw() + xlab("")
# New adjustment which only changes neuronal populations
hseq_ctp_pheno_adjoligo.tmp <- hseq_ctp_pheno.tmp %>%
mutate(sum_neuro_oligo = hseq_Exc + hseq_Inh + hseq_Oligo) %>%
mutate(Exc_Inh_ratio = (hseq_Exc/(hseq_Exc + hseq_Inh))) %>%
mutate(Inh_Exc_ratio = (hseq_Inh/(hseq_Exc + hseq_Inh))) %>%
mutate(hseq_Exc = Exc_Inh_ratio * sum_neuro_oligo) %>%
mutate(hseq_Inh = Inh_Exc_ratio * sum_neuro_oligo) %>%
dplyr::select(-c("hseq_Oligo", "hseq_Neuron", "hseq_Glial"))
# Old adjustment which changed ALL CTPs
# methylcc_ctp_pheno_adjoligo.tmp$sum_nonoligo <- methylcc_ctp_pheno_adjoligo.tmp$Exc + methylcc_ctp_pheno_adjoligo.tmp $Inh + methylcc_ctp_pheno_adjoligo.tmp$NonN_Astro_FGF3R + methylcc_ctp_pheno_adjoligo.tmp$NonN_Micro.Endo_TYROBP
# methylcc_ctp_pheno_adjoligo.tmp$Exc <- methylcc_ctp_pheno_adjoligo.tmp$Exc/methylcc_ctp_pheno_adjoligo.tmp$sum_nonoligo
# methylcc_ctp_pheno_adjoligo.tmp$Inh <- methylcc_ctp_pheno_adjoligo.tmp$Inh/methylcc_ctp_pheno_adjoligo.tmp$sum_nonoligo
# methylcc_ctp_pheno_adjoligo.tmp$NonN_Astro_FGF3R <- methylcc_ctp_pheno_adjoligo.tmp$NonN_Astro_FGF3R/methylcc_ctp_pheno_adjoligo.tmp$sum_nonoligo
# methylcc_ctp_pheno_adjoligo.tmp$NonN_Micro.Endo_TYROBP <- methylcc_ctp_pheno_adjoligo.tmp$NonN_Micro.Endo_TYROBP/methylcc_ctp_pheno_adjoligo.tmp$sum_nonoligo
hseq_ctp_pheno_adjoligo.tmp$hseq_Neuron <- rowSums(hseq_ctp_pheno_adjoligo.tmp[,neuron_n])
hseq_ctp_pheno_adjoligo.tmp$hseq_Glial <- rowSums(hseq_ctp_pheno_adjoligo.tmp[,glial_rmoligo_n])
if (write_out == TRUE) {
write.table(hseq_ctp_pheno_adjoligo.tmp, paste(out_dir, "/hseq_adjoligo_ctp_pheno.raw.txt", sep = ""), col.names = T, row.names = F, quote = F, sep = "\t")
}
hseq_ctp_pheno_adjoligo.long <- melt(hseq_ctp_pheno_adjoligo.tmp, measure.vars = c(all_of(neuron_n), all_of(glial_rmoligo_n), "hseq_Neuron", "hseq_Glial"), variable.name = "celltype", value.name = "CTP")
hseq_ctp_pheno_adjoligo.long$Level <- ifelse(hseq_ctp_pheno_adjoligo.long$celltype %in% c("hseq_Neuron", "hseq_Glial"), "Level1", "Level2")
hseq_ctp_pheno_adjoligo.long$Region <- "PFC"
hseq_ctp_pheno_adjoligo.long$celltype <- gsub("hseq_", "", hseq_ctp_pheno_adjoligo.long$celltype)
hseq_ctp_pheno_adjoligo.long %>% mutate(celltype = fct_reorder(celltype, match(celltype, rev(c("Exc", "Inh", "Astro", "Endo", "Micro", "OPC"))))) %>%
ggplot(aes(x=celltype, y = CTP, fill = group)) +
geom_boxplot(position = position_dodge(.8),alpha=.8,outlier.size = 0)+
geom_point(position = position_jitterdodge(.1),alpha=.3) +
facet_grid(Level~Region, scales='free', space = 'free_y') + coord_flip() + theme_bw() + xlab("")
ctp_cols_mcc <- gsub("hseq", "mcc", ctp_cols)
neuron_n_mcc <- gsub("hseq", "mcc", neuron_n)
glial_rmoligo_n_mcc <- gsub("hseq", "mcc", glial_rmoligo_n)
# - boxplot
methylcc_ctp_pheno.tmp <- ctp_pheno %>% dplyr::select(c(all_of(keep_cols), all_of(ctp_cols_mcc)))
if (write_out == TRUE) {
write.table(methylcc_ctp_pheno.tmp, paste(out_dir, "/methylcc_ctp_pheno.raw.txt", sep = ""), col.names = T, row.names = F, quote = F, sep = "\t")
}
methylcc_ctp_pheno.long <- melt(methylcc_ctp_pheno.tmp, measure.vars = c(all_of(ctp_cols_mcc)), variable.name = "celltype", value.name = "CTP")
methylcc_ctp_pheno.long$Level <- ifelse(methylcc_ctp_pheno.long$celltype %in% c("mcc_Neuron", "mcc_Glial"), "Level1", "Level2")
methylcc_ctp_pheno.long$Region <- "PFC"
methylcc_ctp_pheno.long$celltype <- gsub("mcc_", "", methylcc_ctp_pheno.long$celltype)
methylcc_ctp_pheno.long %>% mutate(celltype = fct_reorder(celltype, match(celltype, rev(c("Exc", "Inh", "Astro", "Endo", "Micro", "Oligo", "OPC"))))) %>%
ggplot(aes(x=celltype, y = CTP, fill = group)) +
geom_boxplot(position = position_dodge(.8),alpha=.8,outlier.size = 0)+
geom_point(position = position_jitterdodge(.1),alpha=.3) +
facet_grid(Level~Region, scales='free', space = 'free_y') + coord_flip() + theme_bw() + xlab("")
# New adjustment which only changes neuronal populations
methylcc_ctp_pheno_adjoligo.tmp <- methylcc_ctp_pheno.tmp %>%
mutate(sum_neuro_oligo = mcc_Exc + mcc_Inh + mcc_Oligo) %>%
mutate(Exc_Inh_ratio = (mcc_Exc/(mcc_Exc + mcc_Inh))) %>%
mutate(Inh_Exc_ratio = (mcc_Inh/(mcc_Exc + mcc_Inh))) %>%
mutate(mcc_Exc = Exc_Inh_ratio * sum_neuro_oligo) %>%
mutate(mcc_Inh = Inh_Exc_ratio * sum_neuro_oligo) %>%
dplyr::select(-c("mcc_Oligo", "mcc_Neuron", "mcc_Glial"))
# Old adjustment which changed ALL CTPs
# methylcc_ctp_pheno_adjoligo.tmp$sum_nonoligo <- methylcc_ctp_pheno_adjoligo.tmp$Exc + methylcc_ctp_pheno_adjoligo.tmp $Inh + methylcc_ctp_pheno_adjoligo.tmp$NonN_Astro_FGF3R + methylcc_ctp_pheno_adjoligo.tmp$NonN_Micro.Endo_TYROBP
# methylcc_ctp_pheno_adjoligo.tmp$Exc <- methylcc_ctp_pheno_adjoligo.tmp$Exc/methylcc_ctp_pheno_adjoligo.tmp$sum_nonoligo
# methylcc_ctp_pheno_adjoligo.tmp$Inh <- methylcc_ctp_pheno_adjoligo.tmp$Inh/methylcc_ctp_pheno_adjoligo.tmp$sum_nonoligo
# methylcc_ctp_pheno_adjoligo.tmp$NonN_Astro_FGF3R <- methylcc_ctp_pheno_adjoligo.tmp$NonN_Astro_FGF3R/methylcc_ctp_pheno_adjoligo.tmp$sum_nonoligo
# methylcc_ctp_pheno_adjoligo.tmp$NonN_Micro.Endo_TYROBP <- methylcc_ctp_pheno_adjoligo.tmp$NonN_Micro.Endo_TYROBP/methylcc_ctp_pheno_adjoligo.tmp$sum_nonoligo
methylcc_ctp_pheno_adjoligo.tmp$mcc_Neuron <- rowSums(methylcc_ctp_pheno_adjoligo.tmp[,neuron_n_mcc])
methylcc_ctp_pheno_adjoligo.tmp$mcc_Glial <- rowSums(methylcc_ctp_pheno_adjoligo.tmp[,glial_rmoligo_n_mcc])
if (write_out == TRUE) {
write.table(methylcc_ctp_pheno_adjoligo.tmp, paste(out_dir, "/methylcc_adjoligo_ctp_pheno.raw.txt", sep = ""), col.names = T, row.names = F, quote = F, sep = "\t")
}
methylcc_ctp_pheno_adjoligo.long <- melt(methylcc_ctp_pheno_adjoligo.tmp, measure.vars = c(all_of(neuron_n_mcc), all_of(glial_rmoligo_n_mcc), "mcc_Neuron", "mcc_Glial"), variable.name = "celltype", value.name = "CTP")
methylcc_ctp_pheno_adjoligo.long$Level <- ifelse(methylcc_ctp_pheno_adjoligo.long$celltype %in% c("mcc_Neuron", "mcc_Glial"), "Level1", "Level2")
methylcc_ctp_pheno_adjoligo.long$Region <- "PFC"
methylcc_ctp_pheno_adjoligo.long$celltype <- gsub("mcc_", "", methylcc_ctp_pheno_adjoligo.long$celltype)
methylcc_ctp_pheno_adjoligo.long %>% mutate(celltype = fct_reorder(celltype, match(celltype, rev(c("Exc", "Inh", "Astro", "Endo", "Micro", "OPC"))))) %>%
ggplot(aes(x=celltype, y = CTP, fill = group)) +
geom_boxplot(position = position_dodge(.8),alpha=.8,outlier.size = 0)+
geom_point(position = position_jitterdodge(.1),alpha=.3) +
facet_grid(Level~Region, scales='free', space = 'free_y') + coord_flip() + theme_bw() + xlab("")
# - boxplot
coefs_sub_ctp_pheno.tmp <- ctp_pheno %>% dplyr::select(c(all_of(keep_cols), "h_GLU", "h_GABA", "h_ASTRO", "h_ENDO", "h_MONO", "h_OPC", "h_OLIGO"))
coefs_sub_ctp_pheno.tmp$neuron_sum <- coefs_sub_ctp_pheno.tmp$h_GABA + coefs_sub_ctp_pheno.tmp$h_GLU
coefs_sub_ctp_pheno.tmp$glia_sum <- coefs_sub_ctp_pheno.tmp$h_ASTRO + coefs_sub_ctp_pheno.tmp$h_OPC + coefs_sub_ctp_pheno.tmp$h_OLIGO + coefs_sub_ctp_pheno.tmp$h_ENDO + coefs_sub_ctp_pheno.tmp$h_MONO
if (write_out == TRUE) {
write.table(coefs_sub_ctp_pheno.tmp, paste(out_dir, "/coefs_sub_ctp_pheno.raw.txt", sep = ""), col.names = T, row.names = F, quote = F, sep = "\t")
}
coefs_sub_ctp_pheno.long <- melt(coefs_sub_ctp_pheno.tmp, measure.vars = c("h_GLU", "h_GABA", "h_ASTRO", "h_ENDO", "h_MONO", "h_OPC", "h_OLIGO", "neuron_sum", "glia_sum"), variable.name = "celltype", value.name = "CTP")
coefs_sub_ctp_pheno.long$Level <- ifelse(coefs_sub_ctp_pheno.long$celltype %in% c("neuron_sum", "glia_sum"), "Level1", "Level2")
coefs_sub_ctp_pheno.long$Region <- "PFC"
ggplot(coefs_sub_ctp_pheno.long, aes(x=celltype, y = CTP, fill = group)) +
geom_boxplot(position = position_dodge(.8),alpha=.8,outlier.size = 0)+
geom_point(position = position_jitterdodge(.1),alpha=.3) +
facet_grid(Level~Region, scales='free', space = 'free_y') + coord_flip() + theme_bw() + xlab("")
keep_pheno_col <- hseq_ctp_pheno.tmp %>% dplyr::select(all_of(keep_cols))
hseq_ctp_pheno.tmp.clr <- as.matrix(hseq_ctp_pheno.tmp %>% dplyr::select(c(all_of(level2_cols))))
# Apply zCompositions to empirically get the best replacement of zeros
# - it looks like you have to transpose the matrix (with samples as columns) to do this properly (?)
length(which(hseq_ctp_pheno.tmp.clr == 0)) # check there are no 0s
## [1] 2
check <- which(hseq_ctp_pheno.tmp.clr == 0, arr.ind = T)
hseq_ctp_pheno.tmp.clr_t <- t(hseq_ctp_pheno.tmp.clr)
if (length(which(hseq_ctp_pheno.tmp.clr == 0)) > 0) {
check <- which(hseq_ctp_pheno.tmp.clr_t == 0, arr.ind = T)
hseq_ctp_pheno.tmp.clr0_t <- cmultRepl(hseq_ctp_pheno.tmp.clr_t, output = "p-counts")
hseq_ctp_pheno.tmp.clr0_t[check]
} else {
print("no zcompositions applied as no 0 values")
hseq_ctp_pheno.tmp.clr0_t <- hseq_ctp_pheno.tmp.clr_t
}
## No. corrected values: 2
## [1] 0.003802408 0.003950411
# Perform clr-transform
hseq_ctp_pheno.tmp.clr0.tr <- clr(t(hseq_ctp_pheno.tmp.clr_t))
hseq_ctp_pheno.clr.zc <- data.frame(keep_pheno_col, hseq_ctp_pheno.tmp.clr0.tr)
if (write_out == TRUE) {
write.table(hseq_ctp_pheno.clr.zc, paste(out_dir, "/hseq_ctp_pheno.clr.zc.txt", sep = ""), col.names = T, row.names = F, quote = F, sep = "\t")
}
Make dataframe
hseq_ctp_pheno.clr.zc.long <- melt(hseq_ctp_pheno.clr.zc, measure.vars = c(all_of(level2_cols)), variable.name = "celltype", value.name = "CTP")
#methylcc_ctp_pheno.clr.zc.long$Level <- ifelse(methylcc_ctp_pheno.clr.zc.long$celltype %in% c("neuron_sum", "glia_sum"), "Level1", "Level2")
hseq_ctp_pheno.clr.zc.long$Region <- "PFC"
ggplot(hseq_ctp_pheno.clr.zc.long, aes(x=celltype, y = CTP, fill = group)) +
geom_boxplot(position = position_dodge(.8),alpha=.8,outlier.size = 0)+
geom_point(position = position_jitterdodge(.1),alpha=.3) +
coord_flip() + theme_bw() + xlab("")
ggtitle("zCompositions")
## $title
## [1] "zCompositions"
##
## attr(,"class")
## [1] "labels"
#facet_grid(Region, scales='free', space = 'free_y') + coord_flip() + theme_bw() + xlab("")
Justification: does it make sense to run an ANOVA on CTP ~ group, if the CTPs are related?
Eg. what does variance analysis look like when the y-variable for each individual is non-independent?
N.B. IDs go as columns, CTPs as rows
keep_pheno_col <- hseq_ctp_pheno.tmp %>% dplyr::select(all_of(keep_cols))
hseq_ctp_pheno.tmp.clr <- as.matrix(hseq_ctp_pheno.tmp %>% dplyr::select(c(all_of(level2_cols))))
# Make minimum CTP = 1e-3
length(which(hseq_ctp_pheno.tmp.clr <= 1e-3)) # check there are no 0s
## [1] 4
hseq_ctp_pheno.tmp.clr0 <- hseq_ctp_pheno.tmp.clr
hseq_ctp_pheno.tmp.clr0[which(hseq_ctp_pheno.tmp.clr0 <= 1e-3)] <- 1e-3
# Perform clr-transform
hseq_ctp_pheno.tmp.clr0.tr <- clr(hseq_ctp_pheno.tmp.clr0)
hseq_ctp_pheno.clr.1e3 <- data.frame(keep_pheno_col, hseq_ctp_pheno.tmp.clr0.tr)
if (write_out == TRUE) {
write.table(hseq_ctp_pheno.clr.1e3, paste(out_dir, "/hseq_ctp_pheno.clr.1e3.txt", sep = ""), col.names = T, row.names = F, quote = F, sep = "\t")
}
Make dataframe
hseq_ctp_pheno.clr.1e3.long <- melt(hseq_ctp_pheno.clr.1e3, measure.vars = c(all_of(level2_cols)), variable.name = "celltype", value.name = "CTP")
hseq_ctp_pheno.clr.1e3.long$Region <- "PFC"
ggplot(hseq_ctp_pheno.clr.1e3.long, aes(x=celltype, y = CTP, fill = group)) +
geom_boxplot(position = position_dodge(.8),alpha=.8,outlier.size = 0)+
geom_point(position = position_jitterdodge(.1),alpha=.3) +
coord_flip() + theme_bw() + xlab("") +
ggtitle("Minimum value = 1e-3")
keep_pheno_col <- hseq_ctp_pheno_adjoligo.tmp %>% dplyr::select(all_of(keep_cols))
hseq_ctp_pheno_adjoligo.tmp.clr <- as.matrix(hseq_ctp_pheno_adjoligo.tmp %>% dplyr::select(c(all_of(neuron_n), all_of(glial_rmoligo_n))))
# Make minimum CTP = 1e-3
length(which(hseq_ctp_pheno_adjoligo.tmp.clr <= 1e-3)) # check there are no 0s
## [1] 4
hseq_ctp_pheno_adjoligo.tmp.clr0 <- hseq_ctp_pheno_adjoligo.tmp.clr
hseq_ctp_pheno_adjoligo.tmp.clr0[which(hseq_ctp_pheno_adjoligo.tmp.clr0 <= 1e-3)] <- 1e-3
# Perform clr-transform
hseq_ctp_pheno_adjoligo.tmp.clr0.tr <- clr(hseq_ctp_pheno_adjoligo.tmp.clr0)
hseq_ctp_pheno_adjoligo.clr.1e3 <- data.frame(keep_pheno_col, hseq_ctp_pheno_adjoligo.tmp.clr0.tr)
if (write_out == TRUE) {
write.table(hseq_ctp_pheno_adjoligo.clr.1e3, paste(out_dir, "/hseq_adjoligo_ctp_pheno.clr.1e3.txt", sep = ""), col.names = T, row.names = F, quote = F, sep = "\t")
}
Make dataframe
hseq_ctp_pheno_adjoligo.clr.1e3.long <- melt(hseq_ctp_pheno_adjoligo.clr.1e3, measure.vars = c(all_of(neuron_n), all_of(glial_rmoligo_n)), variable.name = "celltype", value.name = "CTP")
#methylcc_ctp_pheno.clr.zc.long$Level <- ifelse(methylcc_ctp_pheno.clr.zc.long$celltype %in% c("neuron_sum", "glia_sum"), "Level1", "Level2")
hseq_ctp_pheno_adjoligo.clr.1e3.long$Region <- "PFC"
ggplot(hseq_ctp_pheno_adjoligo.clr.1e3.long, aes(x=celltype, y = CTP, fill = group)) +
geom_boxplot(position = position_dodge(.8),alpha=.8,outlier.size = 0)+
geom_point(position = position_jitterdodge(.1),alpha=.3) +
coord_flip() + theme_bw() + xlab("") +
ggtitle("Minimum value = 1e-3")
#facet_grid(Region, scales='free', space = 'free_y') + coord_flip() + theme_bw() + xlab("")
if (cellnum == 7) {
summary(lm(cbind(hseq_Exc, hseq_Inh, hseq_Astro, hseq_Endo, hseq_Micro, hseq_Oligo, hseq_OPC) ~ group, data = hseq_ctp_pheno.tmp))
summary(lm(cbind(hseq_Exc, hseq_Inh, hseq_Astro, hseq_Endo, hseq_Micro, hseq_Oligo, hseq_OPC) ~ group + age + sex + batch + race + pmi, data = hseq_ctp_pheno.tmp))
} else if (cellnum == 9) {
summary(lm(cbind(Exc_upper, Exc_lower, Inh_MGE, Inh_CGE, NonN_Astro_FGF3R, NonN_Endo, NonN_Micro, NonN_Oligo_MBP, NonN_OPC) ~ group, data = methylcc_ctp_pheno.tmp))
summary(lm(cbind(Exc_upper, Exc_lower, Inh_MGE, Inh_CGE, NonN_Astro_FGF3R, NonN_Endo, NonN_Micro, NonN_Oligo_MBP, NonN_OPC) ~ group + age + sex + batch + race + pmi, data = methylcc_ctp_pheno.tmp))
}
## Response hseq_Exc :
##
## Call:
## lm(formula = hseq_Exc ~ group + age + sex + batch + race + pmi,
## data = hseq_ctp_pheno.tmp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.107495 -0.018705 -0.000705 0.019351 0.103232
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.2038934 0.0233193 8.744 < 2e-16 ***
## groupCTL -0.0076950 0.0024354 -3.160 0.001647 **
## age 0.0004897 0.0002560 1.913 0.056194 .
## sexM -0.0094710 0.0024344 -3.891 0.000109 ***
## batch1 -0.0090001 0.0023925 -3.762 0.000183 ***
## race 0.0036690 0.0041822 0.877 0.380633
## pmi 0.0001779 0.0002006 0.887 0.375441
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03075 on 708 degrees of freedom
## Multiple R-squared: 0.06801, Adjusted R-squared: 0.06011
## F-statistic: 8.611 on 6 and 708 DF, p-value: 4.688e-09
##
##
## Response hseq_Inh :
##
## Call:
## lm(formula = hseq_Inh ~ group + age + sex + batch + race + pmi,
## data = hseq_ctp_pheno.tmp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.066291 -0.011975 0.000034 0.012372 0.056259
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.673e-01 1.366e-02 12.246 <2e-16 ***
## groupCTL -2.709e-03 1.427e-03 -1.899 0.058 .
## age -7.524e-05 1.500e-04 -0.502 0.616
## sexM 2.067e-03 1.426e-03 1.449 0.148
## batch1 1.250e-02 1.402e-03 8.918 <2e-16 ***
## race 3.055e-03 2.450e-03 1.247 0.213
## pmi 1.157e-04 1.175e-04 0.985 0.325
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01801 on 708 degrees of freedom
## Multiple R-squared: 0.1129, Adjusted R-squared: 0.1054
## F-statistic: 15.02 on 6 and 708 DF, p-value: 3.199e-16
##
##
## Response hseq_Astro :
##
## Call:
## lm(formula = hseq_Astro ~ group + age + sex + batch + race +
## pmi, data = hseq_ctp_pheno.tmp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.066687 -0.013257 0.000935 0.013897 0.081743
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.1721458 0.0163100 10.555 < 2e-16 ***
## groupCTL 0.0033298 0.0017033 1.955 0.05099 .
## age -0.0000502 0.0001791 -0.280 0.77933
## sexM 0.0026255 0.0017026 1.542 0.12352
## batch1 -0.0108697 0.0016734 -6.496 1.56e-10 ***
## race 0.0041251 0.0029251 1.410 0.15891
## pmi -0.0004393 0.0001403 -3.131 0.00181 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02151 on 708 degrees of freedom
## Multiple R-squared: 0.08521, Adjusted R-squared: 0.07745
## F-statistic: 10.99 on 6 and 708 DF, p-value: 9.917e-12
##
##
## Response hseq_Endo :
##
## Call:
## lm(formula = hseq_Endo ~ group + age + sex + batch + race + pmi,
## data = hseq_ctp_pheno.tmp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0136082 -0.0029819 -0.0000122 0.0028257 0.0164040
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.226e-02 3.440e-03 9.378 < 2e-16 ***
## groupCTL 1.651e-03 3.593e-04 4.595 5.12e-06 ***
## age -4.307e-05 3.777e-05 -1.140 0.2546
## sexM -8.341e-04 3.591e-04 -2.323 0.0205 *
## batch1 2.264e-03 3.529e-04 6.416 2.56e-10 ***
## race 8.428e-04 6.170e-04 1.366 0.1723
## pmi -4.548e-05 2.959e-05 -1.537 0.1248
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.004536 on 708 degrees of freedom
## Multiple R-squared: 0.09086, Adjusted R-squared: 0.08315
## F-statistic: 11.79 on 6 and 708 DF, p-value: 1.253e-12
##
##
## Response hseq_Micro :
##
## Call:
## lm(formula = hseq_Micro ~ group + age + sex + batch + race +
## pmi, data = hseq_ctp_pheno.tmp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.016100 -0.003552 -0.000580 0.003322 0.034024
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.766e-02 4.379e-03 6.316 4.73e-10 ***
## groupCTL 6.617e-04 4.574e-04 1.447 0.148
## age -5.379e-05 4.809e-05 -1.119 0.264
## sexM 8.449e-04 4.572e-04 1.848 0.065 .
## batch1 2.122e-03 4.493e-04 4.723 2.80e-06 ***
## race 2.054e-04 7.854e-04 0.262 0.794
## pmi -2.887e-05 3.767e-05 -0.766 0.444
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.005775 on 708 degrees of freedom
## Multiple R-squared: 0.0417, Adjusted R-squared: 0.03358
## F-statistic: 5.134 on 6 and 708 DF, p-value: 3.54e-05
##
##
## Response hseq_Oligo :
##
## Call:
## lm(formula = hseq_Oligo ~ group + age + sex + batch + race +
## pmi, data = hseq_ctp_pheno.tmp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.141722 -0.034924 -0.006579 0.026341 0.234953
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3194024 0.0378681 8.435 <2e-16 ***
## groupCTL 0.0047160 0.0039548 1.192 0.2335
## age -0.0002076 0.0004158 -0.499 0.6177
## sexM 0.0051453 0.0039532 1.302 0.1935
## batch1 -0.0038878 0.0038852 -1.001 0.3173
## race -0.0134220 0.0067915 -1.976 0.0485 *
## pmi 0.0002422 0.0003258 0.744 0.4574
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04994 on 708 degrees of freedom
## Multiple R-squared: 0.01342, Adjusted R-squared: 0.00506
## F-statistic: 1.605 on 6 and 708 DF, p-value: 0.1429
##
##
## Response hseq_OPC :
##
## Call:
## lm(formula = hseq_OPC ~ group + age + sex + batch + race + pmi,
## data = hseq_ctp_pheno.tmp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0089936 -0.0021805 -0.0000361 0.0020671 0.0142636
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.948e-03 2.569e-03 2.315 0.0209 *
## groupCTL 2.833e-04 2.683e-04 1.056 0.2914
## age 2.225e-05 2.821e-05 0.789 0.4306
## sexM -6.107e-05 2.682e-04 -0.228 0.8200
## batch1 4.260e-03 2.636e-04 16.161 <2e-16 ***
## race 9.813e-05 4.608e-04 0.213 0.8314
## pmi -4.323e-05 2.210e-05 -1.956 0.0509 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.003388 on 708 degrees of freedom
## Multiple R-squared: 0.2709, Adjusted R-squared: 0.2648
## F-statistic: 43.85 on 6 and 708 DF, p-value: < 2.2e-16
summary(lm(hseq_Glial ~ group, data = ctp_pheno))
##
## Call:
## lm(formula = hseq_Glial ~ group, data = ctp_pheno)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.109356 -0.026825 -0.002075 0.023073 0.179581
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.516203 0.002347 219.939 < 2e-16 ***
## groupCTL 0.011893 0.003077 3.865 0.000121 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04058 on 713 degrees of freedom
## Multiple R-squared: 0.02052, Adjusted R-squared: 0.01915
## F-statistic: 14.94 on 1 and 713 DF, p-value: 0.0001212
summary(lm(mcc_Glial ~ group, data = ctp_pheno))
##
## Call:
## lm(formula = mcc_Glial ~ group, data = ctp_pheno)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.29512 -0.08984 -0.00749 0.08323 0.35525
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.456096 0.007127 63.996 < 2e-16 ***
## groupCTL 0.024421 0.009344 2.614 0.00915 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1232 on 713 degrees of freedom
## Multiple R-squared: 0.00949, Adjusted R-squared: 0.008101
## F-statistic: 6.831 on 1 and 713 DF, p-value: 0.009147
summary(lm(h_NeuN_neg ~ group, data = ctp_pheno))
##
## Call:
## lm(formula = h_NeuN_neg ~ group, data = ctp_pheno)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.157943 -0.035746 -0.000064 0.033708 0.192851
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.427636 0.003085 138.599 < 2e-16 ***
## groupCTL 0.014123 0.004045 3.491 0.00051 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05335 on 713 degrees of freedom
## Multiple R-squared: 0.01681, Adjusted R-squared: 0.01543
## F-statistic: 12.19 on 1 and 713 DF, p-value: 0.0005102
summary(lm(hseq_Glial ~ group + age + sex + batch + race + pmi, data = ctp_pheno))
##
## Call:
## lm(formula = hseq_Glial ~ group + age + sex + batch + race +
## pmi, data = ctp_pheno)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.112344 -0.026521 -0.001446 0.024434 0.191534
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.5574177 0.0305564 18.242 < 2e-16 ***
## groupCTL 0.0106417 0.0031912 3.335 0.000898 ***
## age -0.0003325 0.0003355 -0.991 0.322078
## sexM 0.0077205 0.0031899 2.420 0.015758 *
## batch1 -0.0061108 0.0031350 -1.949 0.051666 .
## race -0.0081505 0.0054802 -1.487 0.137386
## pmi -0.0003146 0.0002629 -1.197 0.231716
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04029 on 708 degrees of freedom
## Multiple R-squared: 0.04126, Adjusted R-squared: 0.03314
## F-statistic: 5.078 on 6 and 708 DF, p-value: 4.076e-05
summary(lm(mcc_Glial ~ group + age + sex + batch + race + pmi, data = ctp_pheno))
##
## Call:
## lm(formula = mcc_Glial ~ group + age + sex + batch + race + pmi,
## data = ctp_pheno)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.31660 -0.08489 0.00118 0.07883 0.37955
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.6040468 0.0923888 6.538 1.19e-10 ***
## groupCTL 0.0181389 0.0096487 1.880 0.0605 .
## age -0.0013827 0.0010144 -1.363 0.1733
## sexM 0.0091320 0.0096447 0.947 0.3440
## batch1 -0.0401530 0.0094789 -4.236 2.58e-05 ***
## race -0.0011935 0.0165695 -0.072 0.9426
## pmi -0.0002321 0.0007948 -0.292 0.7704
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1218 on 708 degrees of freedom
## Multiple R-squared: 0.03878, Adjusted R-squared: 0.03063
## F-statistic: 4.76 on 6 and 708 DF, p-value: 9.075e-05
summary(lm(h_NeuN_neg ~ group + age + sex + batch + race + pmi, data = ctp_pheno))
##
## Call:
## lm(formula = h_NeuN_neg ~ group + age + sex + batch + race +
## pmi, data = ctp_pheno)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.163478 -0.035515 -0.001043 0.033275 0.188151
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.5180070 0.0391117 13.244 < 2e-16 ***
## groupCTL 0.0100405 0.0040847 2.458 0.01421 *
## age -0.0009043 0.0004294 -2.106 0.03557 *
## sexM 0.0131519 0.0040830 3.221 0.00134 **
## batch1 -0.0231003 0.0040128 -5.757 1.28e-08 ***
## race 0.0032781 0.0070145 0.467 0.64041
## pmi -0.0004772 0.0003365 -1.418 0.15657
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05157 on 708 degrees of freedom
## Multiple R-squared: 0.08766, Adjusted R-squared: 0.07993
## F-statistic: 11.34 on 6 and 708 DF, p-value: 4.045e-12
if (cellnum == 7) {
summary(lm(cbind(hseq_Exc, hseq_Inh, hseq_Astro, hseq_Endo, hseq_Micro, hseq_OPC) ~ group, data = hseq_ctp_pheno_adjoligo.tmp))
summary(lm(cbind(hseq_Exc, hseq_Inh, hseq_Astro, hseq_Endo, hseq_Micro, hseq_OPC) ~ group + age + sex + batch + race + pmi, data = hseq_ctp_pheno_adjoligo.tmp))
} else if (cellnum == 9) {
summary(lm(cbind(Exc_upper, Exc_lower, Inh_MGE, Inh_CGE, NonN_Astro_FGF3R, NonN_Endo, NonN_Micro, NonN_OPC) ~ group, data = methylcc_ctp_pheno_adjoligo.tmp))
summary(lm(cbind(Exc_upper, Exc_lower, Inh_MGE, Inh_CGE, NonN_Astro_FGF3R, NonN_Endo, NonN_Micro, NonN_OPC) ~ group + age + sex + batch + race + pmi, data = methylcc_ctp_pheno_adjoligo.tmp))
}
## Response hseq_Exc :
##
## Call:
## lm(formula = hseq_Exc ~ group + age + sex + batch + race + pmi,
## data = hseq_ctp_pheno_adjoligo.tmp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.082089 -0.020570 -0.000001 0.020986 0.128824
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3781132 0.0248053 15.243 < 2e-16 ***
## groupCTL -0.0056948 0.0025906 -2.198 0.0283 *
## age 0.0005690 0.0002724 2.089 0.0371 *
## sexM -0.0101841 0.0025895 -3.933 9.22e-05 ***
## batch1 -0.0189395 0.0025450 -7.442 2.88e-13 ***
## race -0.0040852 0.0044487 -0.918 0.3588
## pmi 0.0002863 0.0002134 1.342 0.1801
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03271 on 708 degrees of freedom
## Multiple R-squared: 0.1099, Adjusted R-squared: 0.1023
## F-statistic: 14.56 on 6 and 708 DF, p-value: 1.024e-15
##
##
## Response hseq_Inh :
##
## Call:
## lm(formula = hseq_Inh ~ group + age + sex + batch + race + pmi,
## data = hseq_ctp_pheno_adjoligo.tmp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.07570 -0.01805 -0.00226 0.01513 0.12163
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.125e-01 2.111e-02 14.800 < 2e-16 ***
## groupCTL 6.328e-06 2.205e-03 0.003 0.997711
## age -3.621e-04 2.318e-04 -1.562 0.118733
## sexM 7.925e-03 2.204e-03 3.596 0.000346 ***
## batch1 1.855e-02 2.166e-03 8.564 < 2e-16 ***
## race -2.613e-03 3.786e-03 -0.690 0.490439
## pmi 2.495e-04 1.816e-04 1.374 0.169882
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02784 on 708 degrees of freedom
## Multiple R-squared: 0.1186, Adjusted R-squared: 0.1111
## F-statistic: 15.88 on 6 and 708 DF, p-value: < 2.2e-16
##
##
## Response hseq_Astro :
##
## Call:
## lm(formula = hseq_Astro ~ group + age + sex + batch + race +
## pmi, data = hseq_ctp_pheno_adjoligo.tmp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.066687 -0.013257 0.000935 0.013897 0.081743
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.1721458 0.0163100 10.555 < 2e-16 ***
## groupCTL 0.0033298 0.0017033 1.955 0.05099 .
## age -0.0000502 0.0001791 -0.280 0.77933
## sexM 0.0026255 0.0017026 1.542 0.12352
## batch1 -0.0108697 0.0016734 -6.496 1.56e-10 ***
## race 0.0041251 0.0029251 1.410 0.15891
## pmi -0.0004393 0.0001403 -3.131 0.00181 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02151 on 708 degrees of freedom
## Multiple R-squared: 0.08521, Adjusted R-squared: 0.07745
## F-statistic: 10.99 on 6 and 708 DF, p-value: 9.917e-12
##
##
## Response hseq_Endo :
##
## Call:
## lm(formula = hseq_Endo ~ group + age + sex + batch + race + pmi,
## data = hseq_ctp_pheno_adjoligo.tmp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0136082 -0.0029819 -0.0000122 0.0028257 0.0164040
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.226e-02 3.440e-03 9.378 < 2e-16 ***
## groupCTL 1.651e-03 3.593e-04 4.595 5.12e-06 ***
## age -4.307e-05 3.777e-05 -1.140 0.2546
## sexM -8.341e-04 3.591e-04 -2.323 0.0205 *
## batch1 2.264e-03 3.529e-04 6.416 2.56e-10 ***
## race 8.428e-04 6.170e-04 1.366 0.1723
## pmi -4.548e-05 2.959e-05 -1.537 0.1248
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.004536 on 708 degrees of freedom
## Multiple R-squared: 0.09086, Adjusted R-squared: 0.08315
## F-statistic: 11.79 on 6 and 708 DF, p-value: 1.253e-12
##
##
## Response hseq_Micro :
##
## Call:
## lm(formula = hseq_Micro ~ group + age + sex + batch + race +
## pmi, data = hseq_ctp_pheno_adjoligo.tmp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.016100 -0.003552 -0.000580 0.003322 0.034024
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.766e-02 4.379e-03 6.316 4.73e-10 ***
## groupCTL 6.617e-04 4.574e-04 1.447 0.148
## age -5.379e-05 4.809e-05 -1.119 0.264
## sexM 8.449e-04 4.572e-04 1.848 0.065 .
## batch1 2.122e-03 4.493e-04 4.723 2.80e-06 ***
## race 2.054e-04 7.854e-04 0.262 0.794
## pmi -2.887e-05 3.767e-05 -0.766 0.444
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.005775 on 708 degrees of freedom
## Multiple R-squared: 0.0417, Adjusted R-squared: 0.03358
## F-statistic: 5.134 on 6 and 708 DF, p-value: 3.54e-05
##
##
## Response hseq_OPC :
##
## Call:
## lm(formula = hseq_OPC ~ group + age + sex + batch + race + pmi,
## data = hseq_ctp_pheno_adjoligo.tmp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0089936 -0.0021805 -0.0000361 0.0020671 0.0142636
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.948e-03 2.569e-03 2.315 0.0209 *
## groupCTL 2.833e-04 2.683e-04 1.056 0.2914
## age 2.225e-05 2.821e-05 0.789 0.4306
## sexM -6.107e-05 2.682e-04 -0.228 0.8200
## batch1 4.260e-03 2.636e-04 16.161 <2e-16 ***
## race 9.813e-05 4.608e-04 0.213 0.8314
## pmi -4.323e-05 2.210e-05 -1.956 0.0509 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.003388 on 708 degrees of freedom
## Multiple R-squared: 0.2709, Adjusted R-squared: 0.2648
## F-statistic: 43.85 on 6 and 708 DF, p-value: < 2.2e-16
if (cellnum == 7) {
summary(lm(cbind(hseq_Exc, hseq_Inh, hseq_Astro, hseq_Endo, hseq_Micro, hseq_Oligo, hseq_OPC) ~ group, data = hseq_ctp_pheno.clr.1e3))
summary(lm(cbind(hseq_Exc, hseq_Inh, hseq_Astro, hseq_Endo, hseq_Micro, hseq_Oligo, hseq_OPC) ~ group + age + sex + batch + race + pmi, data = hseq_ctp_pheno.clr.1e3))
} else if (cellnum == 9) {
summary(lm(cbind(Exc_upper, Exc_lower, Inh_MGE, Inh_CGE, NonN_Astro_FGF3R, NonN_Endo, NonN_Micro, NonN_Oligo_MBP, NonN_OPC) ~ group, data = methylcc_ctp_pheno.clr.1e3))
summary(lm(cbind(Exc_upper, Exc_lower, Inh_MGE, Inh_CGE, NonN_Astro_FGF3R, NonN_Endo, NonN_Micro, NonN_Oligo_MBP, NonN_OPC) ~ group + age + sex + batch + race + pmi, data = methylcc_ctp_pheno.clr.1e3))
}
## Response hseq_Exc :
##
## Call:
## lm(formula = hseq_Exc ~ group + age + sex + batch + race + pmi,
## data = hseq_ctp_pheno.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.61199 -0.11511 -0.00003 0.11305 0.59511
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.065603 0.143838 7.408 3.65e-13 ***
## groupCTL -0.044809 0.015022 -2.983 0.00295 **
## age 0.002048 0.001579 1.297 0.19516
## sexM -0.043744 0.015016 -2.913 0.00369 **
## batch1 -0.131548 0.014758 -8.914 < 2e-16 ***
## race 0.009632 0.025797 0.373 0.70899
## pmi 0.002068 0.001237 1.671 0.09509 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1897 on 708 degrees of freedom
## Multiple R-squared: 0.1259, Adjusted R-squared: 0.1185
## F-statistic: 17 on 6 and 708 DF, p-value: < 2.2e-16
##
##
## Response hseq_Inh :
##
## Call:
## lm(formula = hseq_Inh ~ group + age + sex + batch + race + pmi,
## data = hseq_ctp_pheno.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.35463 -0.06711 -0.01175 0.05923 0.35623
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.8665104 0.0777024 11.152 < 2e-16 ***
## groupCTL -0.0308171 0.0081149 -3.798 0.000159 ***
## age -0.0005377 0.0008532 -0.630 0.528752
## sexM 0.0079539 0.0081116 0.981 0.327147
## batch1 -0.0200145 0.0079721 -2.511 0.012276 *
## race 0.0121114 0.0139356 0.869 0.385088
## pmi 0.0022198 0.0006684 3.321 0.000943 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1025 on 708 degrees of freedom
## Multiple R-squared: 0.04061, Adjusted R-squared: 0.03248
## F-statistic: 4.995 on 6 and 708 DF, p-value: 5.037e-05
##
##
## Response hseq_Astro :
##
## Call:
## lm(formula = hseq_Astro ~ group + age + sex + batch + race +
## pmi, data = hseq_ctp_pheno.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.52551 -0.09358 -0.00108 0.09982 0.50140
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.9058485 0.1179217 7.682 5.23e-14 ***
## groupCTL 0.0057554 0.0123152 0.467 0.640
## age -0.0005046 0.0012948 -0.390 0.697
## sexM 0.0115732 0.0123102 0.940 0.347
## batch1 -0.1590512 0.0120985 -13.146 < 2e-16 ***
## race 0.0166123 0.0211487 0.785 0.432
## pmi -0.0013803 0.0010144 -1.361 0.174
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1555 on 708 degrees of freedom
## Multiple R-squared: 0.2045, Adjusted R-squared: 0.1978
## F-statistic: 30.34 on 6 and 708 DF, p-value: < 2.2e-16
##
##
## Response hseq_Endo :
##
## Call:
## lm(formula = hseq_Endo ~ group + age + sex + batch + race + pmi,
## data = hseq_ctp_pheno.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.32023 -0.05791 0.00431 0.06553 0.39770
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.846e-01 7.196e-02 -10.903 < 2e-16 ***
## groupCTL 3.906e-02 7.515e-03 5.197 2.66e-07 ***
## age -1.462e-03 7.902e-04 -1.851 0.06461 .
## sexM -3.059e-02 7.512e-03 -4.072 5.19e-05 ***
## batch1 -1.941e-02 7.383e-03 -2.629 0.00876 **
## race 1.936e-02 1.291e-02 1.500 0.13403
## pmi -2.563e-05 6.191e-04 -0.041 0.96699
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09489 on 708 degrees of freedom
## Multiple R-squared: 0.08082, Adjusted R-squared: 0.07303
## F-statistic: 10.38 on 6 and 708 DF, p-value: 4.868e-11
##
##
## Response hseq_Micro :
##
## Call:
## lm(formula = hseq_Micro ~ group + age + sex + batch + race +
## pmi, data = hseq_ctp_pheno.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.69719 -0.08813 -0.00355 0.09737 0.63183
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.9416575 0.1192014 -7.900 1.07e-14 ***
## groupCTL 0.0139691 0.0124489 1.122 0.2622
## age -0.0024015 0.0013088 -1.835 0.0670 .
## sexM 0.0321127 0.0124438 2.581 0.0101 *
## batch1 0.0020239 0.0122298 0.165 0.8686
## race -0.0017069 0.0213782 -0.080 0.9364
## pmi 0.0005018 0.0010254 0.489 0.6247
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1572 on 708 degrees of freedom
## Multiple R-squared: 0.02193, Adjusted R-squared: 0.01364
## F-statistic: 2.646 on 6 and 708 DF, p-value: 0.0152
##
##
## Response hseq_Oligo :
##
## Call:
## lm(formula = hseq_Oligo ~ group + age + sex + batch + race +
## pmi, data = hseq_ctp_pheno.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.79621 -0.13513 -0.01171 0.12703 0.80730
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.520e+00 1.568e-01 9.695 < 2e-16 ***
## groupCTL 2.467e-05 1.637e-02 0.002 0.9988
## age -9.140e-04 1.721e-03 -0.531 0.5955
## sexM 1.194e-02 1.636e-02 0.729 0.4660
## batch1 -1.035e-01 1.608e-02 -6.434 2.29e-10 ***
## race -5.163e-02 2.811e-02 -1.837 0.0667 .
## pmi 2.107e-03 1.348e-03 1.563 0.1186
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2067 on 708 degrees of freedom
## Multiple R-squared: 0.06194, Adjusted R-squared: 0.05399
## F-statistic: 7.791 on 6 and 708 DF, p-value: 3.906e-08
##
##
## Response hseq_OPC :
##
## Call:
## lm(formula = hseq_OPC ~ group + age + sex + batch + race + pmi,
## data = hseq_ctp_pheno.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.57065 -0.12977 0.04475 0.19557 0.81437
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.631386 0.251875 -10.447 <2e-16 ***
## groupCTL 0.016821 0.026305 0.639 0.5227
## age 0.003772 0.002766 1.364 0.1730
## sexM 0.010758 0.026294 0.409 0.6826
## batch1 0.431467 0.025842 16.696 <2e-16 ***
## race -0.004378 0.045173 -0.097 0.9228
## pmi -0.005491 0.002167 -2.534 0.0115 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3321 on 708 degrees of freedom
## Multiple R-squared: 0.2863, Adjusted R-squared: 0.2803
## F-statistic: 47.34 on 6 and 708 DF, p-value: < 2.2e-16
if (cellnum == 7) {
summary(lm(cbind(hseq_Exc, hseq_Inh, hseq_Astro, hseq_Endo, hseq_Micro, hseq_OPC) ~ group, data = hseq_ctp_pheno_adjoligo.clr.1e3))
summary(lm(cbind(hseq_Exc, hseq_Inh, hseq_Astro, hseq_Endo, hseq_Micro, hseq_OPC) ~ group + age + sex + batch + race + pmi, data = hseq_ctp_pheno_adjoligo.clr.1e3))
} else if (cellnum == 9) {
summary(lm(cbind(Exc_upper, Exc_lower, Inh_MGE, Inh_CGE, NonN_Astro_FGF3R, NonN_Endo, NonN_Micro, NonN_OPC) ~ group, data = methylcc_ctp_pheno_adjoligo.clr.1e3))
summary(lm(cbind(Exc_upper, Exc_lower, Inh_MGE, Inh_CGE, NonN_Astro_FGF3R, NonN_Endo, NonN_Micro, NonN_OPC) ~ group + age + sex + batch + race + pmi, data = methylcc_ctp_pheno_adjoligo.clr.1e3))
}
## Response hseq_Exc :
##
## Call:
## lm(formula = hseq_Exc ~ group + age + sex + batch + race + pmi,
## data = hseq_ctp_pheno_adjoligo.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.55069 -0.11510 0.00097 0.10484 0.70219
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.729978 0.138015 12.535 <2e-16 ***
## groupCTL -0.033143 0.014414 -2.299 0.0218 *
## age 0.001445 0.001515 0.954 0.3405
## sexM -0.031724 0.014408 -2.202 0.0280 *
## batch1 -0.154971 0.014160 -10.944 <2e-16 ***
## race -0.016924 0.024752 -0.684 0.4944
## pmi 0.002514 0.001187 2.117 0.0346 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.182 on 708 degrees of freedom
## Multiple R-squared: 0.1582, Adjusted R-squared: 0.151
## F-statistic: 22.17 on 6 and 708 DF, p-value: < 2.2e-16
##
##
## Response hseq_Inh :
##
## Call:
## lm(formula = hseq_Inh ~ group + age + sex + batch + race + pmi,
## data = hseq_ctp_pheno_adjoligo.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.51681 -0.07513 -0.00952 0.06449 0.40447
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.5308855 0.0824923 18.558 < 2e-16 ***
## groupCTL -0.0191510 0.0086151 -2.223 0.026534 *
## age -0.0011403 0.0009058 -1.259 0.208462
## sexM 0.0199741 0.0086116 2.319 0.020655 *
## batch1 -0.0434371 0.0084636 -5.132 3.7e-07 ***
## race -0.0144444 0.0147946 -0.976 0.329236
## pmi 0.0026654 0.0007096 3.756 0.000187 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1088 on 708 degrees of freedom
## Multiple R-squared: 0.06434, Adjusted R-squared: 0.05641
## F-statistic: 8.114 on 6 and 708 DF, p-value: 1.696e-08
##
##
## Response hseq_Astro :
##
## Call:
## lm(formula = hseq_Astro ~ group + age + sex + batch + race +
## pmi, data = hseq_ctp_pheno_adjoligo.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.57673 -0.10045 -0.00202 0.10473 0.56930
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.536e-01 1.261e-01 7.559 1.26e-13 ***
## groupCTL -7.155e-05 1.317e-02 -0.005 0.996
## age -4.318e-04 1.385e-03 -0.312 0.755
## sexM 8.547e-03 1.317e-02 0.649 0.517
## batch1 -1.732e-01 1.294e-02 -13.383 < 2e-16 ***
## race 1.698e-02 2.262e-02 0.751 0.453
## pmi -1.076e-03 1.085e-03 -0.992 0.322
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1663 on 708 degrees of freedom
## Multiple R-squared: 0.2073, Adjusted R-squared: 0.2005
## F-statistic: 30.85 on 6 and 708 DF, p-value: < 2.2e-16
##
##
## Response hseq_Endo :
##
## Call:
## lm(formula = hseq_Endo ~ group + age + sex + batch + race + pmi,
## data = hseq_ctp_pheno_adjoligo.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.33734 -0.05729 0.00388 0.06036 0.46129
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.7368667 0.0712418 -10.343 < 2e-16 ***
## groupCTL 0.0332291 0.0074402 4.466 9.27e-06 ***
## age -0.0013897 0.0007822 -1.777 0.0761 .
## sexM -0.0336160 0.0074371 -4.520 7.25e-06 ***
## batch1 -0.0335642 0.0073093 -4.592 5.20e-06 ***
## race 0.0197309 0.0127769 1.544 0.1230
## pmi 0.0002783 0.0006129 0.454 0.6499
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09394 on 708 degrees of freedom
## Multiple R-squared: 0.09314, Adjusted R-squared: 0.08546
## F-statistic: 12.12 on 6 and 708 DF, p-value: 5.389e-13
##
##
## Response hseq_Micro :
##
## Call:
## lm(formula = hseq_Micro ~ group + age + sex + batch + race +
## pmi, data = hseq_ctp_pheno_adjoligo.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.63884 -0.08541 -0.00038 0.09224 0.58907
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.8939247 0.1096891 -8.150 1.65e-15 ***
## groupCTL 0.0081422 0.0114555 0.711 0.4775
## age -0.0023287 0.0012044 -1.933 0.0536 .
## sexM 0.0290867 0.0114508 2.540 0.0113 *
## batch1 -0.0121320 0.0112539 -1.078 0.2814
## race -0.0013368 0.0196722 -0.068 0.9458
## pmi 0.0008058 0.0009436 0.854 0.3934
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1446 on 708 degrees of freedom
## Multiple R-squared: 0.02286, Adjusted R-squared: 0.01458
## F-statistic: 2.761 on 6 and 708 DF, p-value: 0.01167
##
##
## Response hseq_OPC :
##
## Call:
## lm(formula = hseq_OPC ~ group + age + sex + batch + race + pmi,
## data = hseq_ctp_pheno_adjoligo.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.50274 -0.11903 0.04151 0.18724 0.75837
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.583653 0.240254 -10.754 <2e-16 ***
## groupCTL 0.010994 0.025091 0.438 0.6614
## age 0.003845 0.002638 1.458 0.1454
## sexM 0.007732 0.025081 0.308 0.7580
## batch1 0.417311 0.024650 16.930 <2e-16 ***
## race -0.004008 0.043088 -0.093 0.9259
## pmi -0.005187 0.002067 -2.510 0.0123 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3168 on 708 degrees of freedom
## Multiple R-squared: 0.2924, Adjusted R-squared: 0.2864
## F-statistic: 48.76 on 6 and 708 DF, p-value: < 2.2e-16
#if (cellnum == 7) {
hseq_ctp_pheno.clr.1e3$age2 <- hseq_ctp_pheno.clr.1e3$age^2
summary(lm(hseq_Endo ~ group + age + age2 + sex, data = hseq_ctp_pheno.clr.1e3 %>% filter(batch == 0)))
##
## Call:
## lm(formula = hseq_Endo ~ group + age + age2 + sex, data = hseq_ctp_pheno.clr.1e3 %>%
## filter(batch == 0))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.29579 -0.06593 0.00042 0.07388 0.37897
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.300e+00 1.392e+00 -0.934 0.3513
## groupCTL 1.272e-02 1.376e-02 0.924 0.3563
## age 1.334e-02 3.385e-02 0.394 0.6938
## age2 -9.728e-05 2.051e-04 -0.474 0.6357
## sexM -2.890e-02 1.381e-02 -2.093 0.0373 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1042 on 259 degrees of freedom
## Multiple R-squared: 0.0328, Adjusted R-squared: 0.01787
## F-statistic: 2.196 on 4 and 259 DF, p-value: 0.06985
summary(lm(hseq_Endo ~ group + age + age2 + sex, data = hseq_ctp_pheno.clr.1e3 %>% filter(batch == 1)))
##
## Call:
## lm(formula = hseq_Endo ~ group + age + age2 + sex, data = hseq_ctp_pheno.clr.1e3 %>%
## filter(batch == 1))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.309532 -0.050033 0.004614 0.059244 0.234283
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.1235107 1.1634108 -0.106 0.915501
## groupCTL 0.0545639 0.0087948 6.204 1.26e-09 ***
## age -0.0183923 0.0281175 -0.654 0.513370
## age2 0.0001061 0.0001692 0.627 0.531092
## sexM -0.0328742 0.0088943 -3.696 0.000246 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08875 on 446 degrees of freedom
## Multiple R-squared: 0.1103, Adjusted R-squared: 0.1023
## F-statistic: 13.82 on 4 and 446 DF, p-value: 1.232e-10
summary(lm(hseq_Exc ~ group + age + age2 + sex, data = hseq_ctp_pheno.clr.1e3 %>% filter(batch == 0)))
##
## Call:
## lm(formula = hseq_Exc ~ group + age + age2 + sex, data = hseq_ctp_pheno.clr.1e3 %>%
## filter(batch == 0))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.62405 -0.12198 0.00683 0.13444 0.61025
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.2078059 2.8329846 -1.838 0.0672 .
## groupCTL -0.0393019 0.0280110 -1.403 0.1618
## age 0.1554426 0.0688946 2.256 0.0249 *
## age2 -0.0009298 0.0004174 -2.228 0.0268 *
## sexM -0.0372614 0.0280999 -1.326 0.1860
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.212 on 259 degrees of freedom
## Multiple R-squared: 0.03766, Adjusted R-squared: 0.0228
## F-statistic: 2.534 on 4 and 259 DF, p-value: 0.04076
summary(lm(hseq_Exc ~ group + age + age2 + sex, data = hseq_ctp_pheno.clr.1e3 %>% filter(batch == 1)))
##
## Call:
## lm(formula = hseq_Exc ~ group + age + age2 + sex, data = hseq_ctp_pheno.clr.1e3 %>%
## filter(batch == 1))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.45746 -0.11957 -0.00239 0.10316 0.56615
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.9022386 2.2956591 1.264 0.2068
## groupCTL -0.0460820 0.0173540 -2.655 0.0082 **
## age -0.0450890 0.0554818 -0.813 0.4168
## age2 0.0002846 0.0003339 0.853 0.3944
## sexM -0.0441951 0.0175503 -2.518 0.0121 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1751 on 446 degrees of freedom
## Multiple R-squared: 0.04693, Adjusted R-squared: 0.03838
## F-statistic: 5.49 on 4 and 446 DF, p-value: 0.0002538
summary(lm(hseq_Inh ~ group + age + age2 + sex, data = hseq_ctp_pheno.clr.1e3 %>% filter(batch == 0)))
##
## Call:
## lm(formula = hseq_Inh ~ group + age + age2 + sex, data = hseq_ctp_pheno.clr.1e3 %>%
## filter(batch == 0))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.35142 -0.07523 -0.01818 0.07286 0.35341
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.1403314 1.5425474 -0.739 0.4604
## groupCTL -0.0310897 0.0152519 -2.038 0.0425 *
## age 0.0492412 0.0375128 1.313 0.1905
## age2 -0.0003022 0.0002273 -1.330 0.1848
## sexM -0.0002573 0.0153003 -0.017 0.9866
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1155 on 259 degrees of freedom
## Multiple R-squared: 0.02101, Adjusted R-squared: 0.005893
## F-statistic: 1.39 on 4 and 259 DF, p-value: 0.2379
summary(lm(hseq_Inh ~ group + age + age2 + sex, data = hseq_ctp_pheno.clr.1e3 %>% filter(batch == 1)))
##
## Call:
## lm(formula = hseq_Inh ~ group + age + age2 + sex, data = hseq_ctp_pheno.clr.1e3 %>%
## filter(batch == 1))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.28285 -0.06655 -0.00748 0.06341 0.32875
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.2747229 1.2476600 2.625 0.00897 **
## groupCTL -0.0276583 0.0094316 -2.932 0.00354 **
## age -0.0587926 0.0301536 -1.950 0.05183 .
## age2 0.0003512 0.0001814 1.936 0.05353 .
## sexM 0.0158703 0.0095383 1.664 0.09685 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09517 on 446 degrees of freedom
## Multiple R-squared: 0.03171, Adjusted R-squared: 0.02302
## F-statistic: 3.651 on 4 and 446 DF, p-value: 0.006117
#}
hseq_ctp_pheno.tmp.clr0_ratio <- data.frame(hseq_ctp_pheno.tmp.clr0)
hseq_ctp_pheno.tmp.clr0_ratio$exc_inh <- hseq_ctp_pheno.tmp.clr0_ratio$hseq_Exc/hseq_ctp_pheno.tmp.clr0_ratio$hseq_Inh
hseq_ctp_pheno.tmp.clr0_ratio.df <- data.frame(keep_pheno_col, hseq_ctp_pheno.tmp.clr0_ratio)
summary(lm(exc_inh ~ group, data = hseq_ctp_pheno.tmp.clr0_ratio.df))
summary(lm(exc_inh ~ group + age + sex + batch + race + pmi, data = hseq_ctp_pheno.tmp.clr0_ratio.df))
# Remove NA
foo <- hseq_ctp_pheno.tmp %>% filter(!is.na(age)) %>% filter(!is.na(sex)) %>% filter(!is.na(batch)) %>% filter(!is.na(race)) %>% filter(!is.na(pmi))
conds <- foo$group
aldex_in <- sapply(data.frame(t(foo[,c(all_of(level2_cols))])*10000), as.integer)
rownames(aldex_in) <- level2_cols
colnames(aldex_in) <- foo$IID
x <- aldex.clr(aldex_in, conds, mc.samples=128, denom="all", verbose=F)
## operating in serial mode
## computing center with all features
x.tt <- aldex.ttest(x, paired.test = F)
x.effect <- aldex.effect(x, verbose = F)
x.all <- data.frame(x.tt, x.effect)
tmp <- rownames(x.all)
x.all <- as.data.frame(sapply(x.all, function(x) signif(x, 2)))
rownames(x.all) <- tmp
datatable(x.all)
par(mfrow=c(1,2))
aldex.plot(x.all, type="MA", test="welch", xlab="Log-ratio abundance",ylab="Difference")
aldex.plot(x.all, type="MW", test="welch", xlab="Dispersion",ylab="Difference")
conds <- foo$group
cov.aldex <- foo %>% dplyr::select(c("group", "age", "sex", "batch", "race", "pmi"))
mm <- model.matrix(~., cov.aldex)
x <- aldex.clr(aldex_in, mm, mc.samples=128, denom="all", verbose=F)
x.glm <- aldex.glm(x, mm)
## Warning in if (verbose) message("running tests for each MC instance:"): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## |--
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -(25%)
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -(50%)
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -(75%)
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## Warning in if (verbose == TRUE) numTicks <- progress(i, k, numTicks): the
## condition has length > 1 and only the first element will be used
## -|
x.all2 <- data.frame(x.glm)
tmp <- rownames(x.all2)
x.all2 <- as.data.frame(sapply(x.all2, function(x) signif(x, 2)))
rownames(x.all2) <- tmp
datatable(x.all2)
# Prepare data
metadata <- hseq_ctp_pheno.tmp %>% dplyr::select(c("IID", "group", "age", "sex", "batch", "race", "pmi"))
rownames(metadata) <- metadata$IID
metadata$group <- as.factor(metadata$group)
ancom_features <- t(hseq_ctp_pheno.tmp[,c(all_of(level2_cols))])*10000
colnames(ancom_features) <- hseq_ctp_pheno.tmp$IID
ancom_features <- ancom_features[,which(colnames(ancom_features) %in% metadata$IID)]
# Check variables in the correct order
identical(colnames(ancom_features), metadata$IID)
## [1] TRUE
ancom_table <- feature_table_pre_process(ancom_features, metadata, "IID", group_var = NULL, out_cut = 0.05, zero_cut = 0.90, lib_cut = 10, neg_lb)
# Run tests with covariates (age + sex + diet)
ancom_out_nocov <- ANCOM(ancom_table$feature_table, ancom_table$meta_data, ancom_table$structure_zeros, main_var = "group", p_adj_method = "BH")
ancom_out_nocov$out <- ancom_out_nocov$out[order(ancom_out_nocov$out$W, decreasing = T),]
datatable(ancom_out_nocov$out)
ancom_out_nocov$fig
# Run tests with covariates (age + sex + diet)
ancom_out_cov <- ANCOM(ancom_table$feature_table, ancom_table$meta_data, ancom_table$structure_zeros, main_var = "group", p_adj_method = "BH", adj_formula = "age + sex + batch + race + pmi")
ancom_out_cov$out <- ancom_out_cov$out[order(ancom_out_cov$out$W, decreasing = T),]
datatable(ancom_out_cov$out)
ancom_out_cov$fig
cov_var <- c("IID", "group", "age", "sex", "batch", "Sample_Plate", "race", "pmi")
# Select columns
neunn <- ctp_pheno %>% dplyr::select(c(all_of(cov_var), "h_NeuN_neg", "hseq_Glial", "mcc_Glial", "sSV1", "sSV2", "h_Glial", "hseq_Oligo", "hseq_OPC", "mcc_Oligo", "mcc_OPC", "h_OLIGO", "h_OPC", "celfie_Glial", "celfie_Oligo", "celfie_OPC", "VAEe3", "VAEe9"))
# y = h_NeuN_neg
neunn.long <- melt(neunn, id.var = c(all_of(cov_var), "h_NeuN_neg"), measure.var = c("hseq_Glial", "mcc_Glial", "h_Glial", "celfie_Glial", "sSV2", "VAEe3"), variable.name = c("method"), value.name = "CTP")
ggplot(neunn.long, aes(x=CTP, y=h_NeuN_neg)) + geom_point(aes(colour = Sample_Plate)) + scale_colour_viridis(discrete = TRUE) + geom_smooth(method = "lm") + facet_wrap(~method, scales = "free_x")
## `geom_smooth()` using formula 'y ~ x'
# y: hseq_Glial
neunn.long <- melt(neunn, id.var = c(all_of(cov_var), "hseq_Glial"), measure.var = c("h_NeuN_neg", "mcc_Glial", "h_Glial", "celfie_Glial", "sSV2", "VAEe3"), variable.name = c("method"), value.name = "CTP")
ggplot(neunn.long, aes(x=CTP, y=hseq_Glial)) + geom_point(aes(colour = Sample_Plate)) + scale_colour_viridis(discrete = TRUE) + geom_smooth(method = "lm") + facet_wrap(~method, scales = "free_x")
## `geom_smooth()` using formula 'y ~ x'
# y: mcc_Glial
neunn.long <- melt(neunn, id.var = c(all_of(cov_var), "mcc_Glial"), measure.var = c("h_NeuN_neg", "hseq_Glial", "h_Glial", "celfie_Glial", "sSV2", "VAEe3"), variable.name = c("method"), value.name = "CTP")
ggplot(neunn.long, aes(x=CTP, y=mcc_Glial)) + geom_point(aes(colour = Sample_Plate)) + scale_colour_viridis(discrete = TRUE) + geom_smooth(method = "lm") + facet_wrap(~method, scales = "free_x")
## `geom_smooth()` using formula 'y ~ x'
# y: h_Glial
neunn.long <- melt(neunn, id.var = c(all_of(cov_var), "h_Glial"), measure.var = c("h_NeuN_neg", "hseq_Glial", "mcc_Glial", "celfie_Glial", "sSV2", "VAEe3"), variable.name = c("method"), value.name = "CTP")
ggplot(neunn.long, aes(x=CTP, y=h_Glial)) + geom_point(aes(colour = Sample_Plate)) + scale_colour_viridis(discrete = TRUE) + geom_smooth(method = "lm") + facet_wrap(~method, scales = "free_x")
## `geom_smooth()` using formula 'y ~ x'
# y: celfie_Glial
neunn.long <- melt(neunn, id.var = c(all_of(cov_var), "celfie_Glial"), measure.var = c("h_NeuN_neg", "hseq_Glial", "mcc_Glial", "h_Glial", "sSV2", "VAEe3"), variable.name = c("method"), value.name = "CTP")
ggplot(neunn.long, aes(x=CTP, y=celfie_Glial)) + geom_point(aes(colour = Sample_Plate)) + scale_colour_viridis(discrete = TRUE) + geom_smooth(method = "lm") + facet_wrap(~method, scales = "free_x")
## `geom_smooth()` using formula 'y ~ x'
# y: sSV2
neunn.long <- melt(neunn, id.var = c(all_of(cov_var), "sSV2"), measure.var = c("h_NeuN_neg", "hseq_Oligo", "mcc_Oligo", "h_OLIGO", "celfie_Oligo", "VAEe9"), variable.name = c("method"), value.name = "CTP")
ggplot(neunn.long, aes(x=CTP, y=sSV2)) + geom_point(aes(colour = Sample_Plate)) + scale_colour_viridis(discrete = TRUE) + geom_smooth(method = "lm") + facet_wrap(~method, scales = "free_x")
## `geom_smooth()` using formula 'y ~ x'
# y: sSV2
neunn.long <- melt(neunn, id.var = c(all_of(cov_var), "sSV2"), measure.var = c("h_NeuN_neg", "hseq_Glial", "mcc_Glial", "h_Glial", "celfie_Glial", "hseq_Oligo", "hseq_OPC", "mcc_Oligo", "mcc_OPC", "h_OLIGO", "h_OPC", "celfie_Oligo", "celfie_OPC", "VAEe9"), variable.name = c("method"), value.name = "CTP")
ggplot(neunn.long, aes(x=CTP, y=sSV2)) + geom_point(aes(colour = Sample_Plate)) + scale_colour_viridis(discrete = TRUE) + geom_smooth(method = "lm") + facet_wrap(~method, scales = "free_x")
## `geom_smooth()` using formula 'y ~ x'
# y: VAEe3
neunn.long <- melt(neunn, id.var = c(all_of(cov_var), "VAEe9"), measure.var = c("h_NeuN_neg", "hseq_Oligo", "mcc_Oligo", "h_OLIGO", "celfie_Oligo", "sSV2"), variable.name = c("method"), value.name = "CTP")
ggplot(neunn.long, aes(x=CTP, y=VAEe9)) + geom_point(aes(colour = Sample_Plate)) + scale_colour_viridis(discrete = TRUE) + geom_smooth(method = "lm") + facet_wrap(~method, scales = "free_x")
## `geom_smooth()` using formula 'y ~ x'
# Comparison between houseman extremes vs eBayes
ctp_pheno_sub_glial <- data.frame(ctp_pheno[,c("hseq_Glial", "hseq_Astro", "hseq_Endo", "hseq_Micro", "hseq_Oligo", "hseq_OPC", "h_NeuN_neg", "h_Glial", "h_ASTRO", "h_ENDO", "h_MONO", "h_OLIGO", "h_OPC")])
ggpairs(ctp_pheno_sub_glial)
# Comparison between houseman extremes vs methylCC extremes
ctp_pheno_sub_glial <- data.frame(ctp_pheno[,c("hseq_Glial", "hseq_Astro", "hseq_Endo", "hseq_Micro", "hseq_Oligo", "hseq_OPC", "h_NeuN_neg", "h_Glial", "h_ASTRO", "h_ENDO", "h_MONO", "h_OLIGO", "h_OPC")])
ggpairs(ctp_pheno_sub_glial)
# Comparison between houseman extremes vs CelFIE
ctp_pheno_sub_glial <- data.frame(ctp_pheno[,c("hseq_Glial", "hseq_Astro", "hseq_Endo", "hseq_Micro", "hseq_Oligo", "hseq_OPC", "h_NeuN_neg", "celfie_Glial", "celfie_Astro", "celfie_Endo", "celfie_Micro", "celfie_Oligo", "celfie_OPC")])
ggpairs(ctp_pheno_sub_glial)
# Comparison between houseman array vs methylCC extremes
ctp_pheno_sub_glial <- data.frame(ctp_pheno[,c("h_Glial", "h_ASTRO", "h_ENDO", "h_MONO", "h_OLIGO", "h_OPC", "h_NeuN_neg", "mcc_Glial", "mcc_Astro", "mcc_Endo", "mcc_Micro", "mcc_Oligo", "mcc_OPC")])
ggpairs(ctp_pheno_sub_glial)
# Comparison between houseman array vs methylCC extremes
ctp_pheno_sub_glial <- data.frame(ctp_pheno[,c("h_Glial", "h_ASTRO", "h_ENDO", "h_MONO", "h_OLIGO", "h_OPC", "h_NeuN_neg", "mcc_Glial", "mcc_Astro", "mcc_Endo", "mcc_Micro", "mcc_Oligo", "mcc_OPC")])
ggpairs(ctp_pheno_sub_glial)
# Select columns
neunp <- ctp_pheno %>% dplyr::select(c(all_of(cov_var), "h_NeuN_pos", "hseq_Neuron", "mcc_Neuron", "h_Neuron", "celfie_Neuron", "sSV2", "VAEe3"))
# y: h_NeuN_pos
neunp.long <- melt(neunp, id.var = c(all_of(cov_var), "h_NeuN_pos"), measure.var = c("hseq_Neuron", "mcc_Neuron", "h_Neuron", "celfie_Neuron", "sSV2", "VAEe3"), variable.name = c("method"), value.name = "CTP")
ggplot(neunp.long, aes(x=CTP, y=h_NeuN_pos)) + geom_point(aes(colour = Sample_Plate)) + scale_colour_viridis(discrete = T) + geom_smooth(method = "lm") + facet_wrap(~method, scales = "free_x")
## `geom_smooth()` using formula 'y ~ x'
# y: hseq_Neuron
neunp.long <- melt(neunp, id.var = c(all_of(cov_var), "hseq_Neuron"), measure.var = c("h_NeuN_pos", "mcc_Neuron", "h_Neuron", "celfie_Neuron", "sSV2", "VAEe3"), variable.name = c("method"), value.name = "CTP")
ggplot(neunp.long, aes(x=CTP, y=hseq_Neuron)) + geom_point(aes(colour = Sample_Plate)) + scale_colour_viridis(discrete = T) + geom_smooth(method = "lm") + facet_wrap(~method, scales = "free_x")
## `geom_smooth()` using formula 'y ~ x'
# y: mcc_Neuron
neunp.long <- melt(neunp, id.var = c(all_of(cov_var), "mcc_Neuron"), measure.var = c("h_NeuN_pos", "hseq_Neuron", "h_Neuron", "celfie_Neuron", "sSV2", "VAEe3"), variable.name = c("method"), value.name = "CTP")
ggplot(neunp.long, aes(x=CTP, y=mcc_Neuron)) + geom_point(aes(colour = Sample_Plate)) + scale_colour_viridis(discrete = T) + geom_smooth(method = "lm") + facet_wrap(~method, scales = "free_x")
## `geom_smooth()` using formula 'y ~ x'
# y: h_Neuron
neunp.long <- melt(neunp, id.var = c(all_of(cov_var), "h_Neuron"), measure.var = c("h_NeuN_pos", "hseq_Neuron", "mcc_Neuron", "celfie_Neuron", "sSV2", "VAEe3"), variable.name = c("method"), value.name = "CTP")
ggplot(neunp.long, aes(x=CTP, y=h_Neuron)) + geom_point(aes(colour = Sample_Plate)) + scale_colour_viridis(discrete = T) + geom_smooth(method = "lm") + facet_wrap(~method, scales = "free_x")
## `geom_smooth()` using formula 'y ~ x'
# y: celfie_Neuron
neunp.long <- melt(neunp, id.var = c(all_of(cov_var), "celfie_Neuron"), measure.var = c("h_NeuN_pos", "hseq_Neuron", "mcc_Neuron", "h_Neuron", "sSV2", "VAEe3"), variable.name = c("method"), value.name = "CTP")
ggplot(neunp.long, aes(x=CTP, y=celfie_Neuron)) + geom_point(aes(colour = Sample_Plate)) + scale_colour_viridis(discrete = T) + geom_smooth(method = "lm") + facet_wrap(~method, scales = "free_x")
## `geom_smooth()` using formula 'y ~ x'
# y: sSV2
neunp.long <- melt(neunp, id.var = c(all_of(cov_var), "sSV2"), measure.var = c("h_NeuN_pos", "hseq_Neuron", "mcc_Neuron", "h_Neuron", "celfie_Neuron", "VAEe3"), variable.name = c("method"), value.name = "CTP")
ggplot(neunp.long, aes(x=CTP, y=sSV2)) + geom_point(aes(colour = Sample_Plate)) + scale_colour_viridis(discrete = T) + geom_smooth(method = "lm") + facet_wrap(~method, scales = "free_x")
## `geom_smooth()` using formula 'y ~ x'
# y: VAEe3
neunp.long <- melt(neunp, id.var = c(all_of(cov_var), "VAEe3"), measure.var = c("h_NeuN_pos", "hseq_Neuron", "mcc_Neuron", "h_Neuron", "celfie_Neuron", "sSV2"), variable.name = c("method"), value.name = "CTP")
ggplot(neunp.long, aes(x=CTP, y=VAEe3)) + geom_point(aes(colour = Sample_Plate)) + scale_colour_viridis(discrete = T) + geom_smooth(method = "lm") + facet_wrap(~method, scales = "free_x")
## `geom_smooth()` using formula 'y ~ x'
if (cellnum == 7) {
ctp_pheno_sub_neuron <- data.frame(ctp_pheno[,c("hseq_Neuron", "hseq_Exc", "hseq_Inh", "h_Neuron", "h_GLU", "h_GABA", "mcc_Neuron", "mcc_Exc", "mcc_Inh", "h_NeuN_pos")])
} else if (cellnum == 9) {
ctp_pheno_sub_neuron <- data.frame(ctp_pheno[,c("Neuron", "Exc_upper", "Exc_lower", "Inh_MGE", "Inh_CGE", "h_NeuN_pos", "h_Neuron", "h_GABA", "h_GLU")])
}
ggpairs(ctp_pheno_sub_neuron)
ctp_pheno_sub <- data.frame(ctp_pheno[,c("hseq_Exc", "hseq_Inh", "hseq_Astro", "hseq_Endo", "hseq_Micro", "hseq_Oligo", "hseq_OPC", "h_GLU", "h_GABA", "h_ASTRO", "h_ENDO", "h_MONO", "h_OLIGO", "h_OPC")])
#svg(paste(out_dir, "/", filen, "_ctp_hseq_harray_cor_matrix.svg", sep = ""), height = 6, width = 6)
jpeg(paste(out_dir, "/", filen, "_ctp_hseq_harray_cor_matrix.jpeg", sep = ""), height = 600, width = 600, quality = 100)
g <- ggduo(ctp_pheno_sub, colnames(ctp_pheno_sub)[1:7], colnames(ctp_pheno_sub)[8:14], xlab = "Houseman + sequencing reference", ylab = "Houseman + array reference", title = "ROSMAP")
print(g)
## plot: [1,1] [>-------------------------------------------------] 2% est: 0s
## plot: [1,2] [=>------------------------------------------------] 4% est: 1s
## plot: [1,3] [==>-----------------------------------------------] 6% est: 3s
## plot: [1,4] [===>----------------------------------------------] 8% est: 3s
## plot: [1,5] [====>---------------------------------------------] 10% est: 2s
## plot: [1,6] [=====>--------------------------------------------] 12% est: 2s
## plot: [1,7] [======>-------------------------------------------] 14% est: 2s
## plot: [2,1] [=======>------------------------------------------] 16% est: 2s
## plot: [2,2] [========>-----------------------------------------] 18% est: 2s
## plot: [2,3] [=========>----------------------------------------] 20% est: 2s
## plot: [2,4] [==========>---------------------------------------] 22% est: 2s
## plot: [2,5] [===========>--------------------------------------] 24% est: 2s
## plot: [2,6] [============>-------------------------------------] 27% est: 2s
## plot: [2,7] [=============>------------------------------------] 29% est: 2s
## plot: [3,1] [==============>-----------------------------------] 31% est: 2s
## plot: [3,2] [===============>----------------------------------] 33% est: 2s
## plot: [3,3] [================>---------------------------------] 35% est: 2s
## plot: [3,4] [=================>--------------------------------] 37% est: 1s
## plot: [3,5] [==================>-------------------------------] 39% est: 1s
## plot: [3,6] [===================>------------------------------] 41% est: 1s
## plot: [3,7] [====================>-----------------------------] 43% est: 1s
## plot: [4,1] [=====================>----------------------------] 45% est: 1s
## plot: [4,2] [======================>---------------------------] 47% est: 1s
## plot: [4,3] [=======================>--------------------------] 49% est: 1s
## plot: [4,4] [=========================>------------------------] 51% est: 1s
## plot: [4,5] [==========================>-----------------------] 53% est: 1s
## plot: [4,6] [===========================>----------------------] 55% est: 1s
## plot: [4,7] [============================>---------------------] 57% est: 1s
## plot: [5,1] [=============================>--------------------] 59% est: 1s
## plot: [5,2] [==============================>-------------------] 61% est: 1s
## plot: [5,3] [===============================>------------------] 63% est: 1s
## plot: [5,4] [================================>-----------------] 65% est: 1s
## plot: [5,5] [=================================>----------------] 67% est: 1s
## plot: [5,6] [==================================>---------------] 69% est: 1s
## plot: [5,7] [===================================>--------------] 71% est: 1s
## plot: [6,1] [====================================>-------------] 73% est: 1s
## plot: [6,2] [=====================================>------------] 76% est: 1s
## plot: [6,3] [======================================>-----------] 78% est: 1s
## plot: [6,4] [=======================================>----------] 80% est: 0s
## plot: [6,5] [========================================>---------] 82% est: 0s
## plot: [6,6] [=========================================>--------] 84% est: 0s
## plot: [6,7] [==========================================>-------] 86% est: 0s
## plot: [7,1] [===========================================>------] 88% est: 0s
## plot: [7,2] [============================================>-----] 90% est: 0s
## plot: [7,3] [=============================================>----] 92% est: 0s
## plot: [7,4] [==============================================>---] 94% est: 0s
## plot: [7,5] [===============================================>--] 96% est: 0s
## plot: [7,6] [================================================>-] 98% est: 0s
## plot: [7,7] [==================================================]100% est: 0s
dev.off()
## X11cairo
## 2
if (cellnum == 7) {
ctp_pheno_sv_neuron <- data.frame(ctp_pheno[,c("h_NeuN_pos", "hseq_Neuron", "hseq_Exc", "hseq_Inh", "sSV1", "sSV2", "sSV3", "sSV4", "sSV5", "sSV6", "sSV7", "sSV8")])
} else if (cellnum == 9) {
ctp_pheno_sv_neuron <- data.frame(ctp_pheno[,c("h_NeuN_pos", "Neuron", "Exc_upper", "Exc_lower", "Inh_MGE", "Inh_CGE", "sSV1", "sSV2", "sSV3", "sSV4", "sSV5", "sSV6", "sSV7", "sSV8")])
}
ggpairs(ctp_pheno_sv_neuron)
ctp_pheno_sv_glial <- data.frame(ctp_pheno[,c("h_NeuN_neg", "hseq_Glial", "hseq_Astro", "hseq_Endo", "hseq_Micro", "hseq_Oligo", "hseq_OPC", "sSV1", "sSV2", "sSV3", "sSV4", "sSV5", "sSV6", "sSV7", "sSV8")])
ggpairs(ctp_pheno_sv_glial)
ihc_dir <- "~/shared-gandalm/GenomicDatasets/ROSMAP/ROSMAP_IHC"
ihc_neuro <- fread(paste(ihc_dir, "/IHC.neuro.txt", sep = ""), header = T)
ihc_astro <- fread(paste(ihc_dir, "/IHC.astro.txt", sep = ""), header = T)
ihc_endo <- fread(paste(ihc_dir, "/IHC.endo.txt", sep = ""), header = T)
ihc_micro <- fread(paste(ihc_dir, "/IHC.microglia.txt", sep = ""), header = T)
ihc_oligo <- fread(paste(ihc_dir, "/IHC.oligo.txt", sep = ""), header = T)
ihc <- data.frame(projid = colnames(ihc_neuro), neuro_ihc = t(ihc_neuro[1,])) %>%
inner_join(data.frame(projid = colnames(ihc_astro), astro_ihc = t(ihc_astro[1,])), by = "projid") %>%
inner_join(data.frame(projid = colnames(ihc_endo), endo_ihc = t(ihc_endo[1,])), by = "projid") %>%
inner_join(data.frame(projid = colnames(ihc_micro), micro_ihc = t(ihc_micro[1,])), by = "projid") %>%
inner_join(data.frame(projid = colnames(ihc_oligo), oligo_ihc = t(ihc_oligo[1,])), by = "projid")
ihc$projid <- as.integer(ihc$projid)
ihc <- as.matrix(ihc)
ihc[which(is.na(ihc))] <- 0
ihc <- data.frame(ihc)
ctp_pheno_ihc <- inner_join(ctp_pheno, ihc, by = "projid")
ctp_pheno_ihc_hseq <- ctp_pheno_ihc[,c("IID", level2_cols, "neuro_ihc", "astro_ihc", "endo_ihc", "micro_ihc", "oligo_ihc")]
# Make neuron aggregated variable for comparison to IHC
# - "transfer" all oligodendrocyte fraction to the Neuron fraction
ctp_pheno_ihc_hseq_rmoligo <- ctp_pheno_ihc_hseq %>%
mutate(hseq_Neuron = hseq_Exc + hseq_Inh + hseq_Oligo) %>%
mutate(neuro_ihc = neuro_ihc + oligo_ihc) %>%
dplyr::select(-c("hseq_Oligo", "oligo_ihc"))
# clr transform
ihc.clr.1e3.tmp0 <- as.matrix(ihc[,c("neuro_ihc", "astro_ihc", "endo_ihc", "micro_ihc", "oligo_ihc")])
ihc.clr.1e3.tmp0[which(ihc.clr.1e3.tmp0 < 1e-3)] <- 1e-3
ihc.clr.1e3.tmp <- clr(ihc.clr.1e3.tmp0)
ihc.clr.1e3 <- data.frame(projid = ihc$projid, ihc.clr.1e3.tmp)
# Redo the hseq clr-transform with Neuron sum
hseq_ctp_pheno.tmp.clr <- as.matrix(hseq_ctp_pheno.tmp %>% dplyr::select(c("hseq_Neuron", "hseq_Astro", "hseq_Endo", "hseq_Micro", "hseq_Oligo", "hseq_OPC")))
length(which(hseq_ctp_pheno.tmp.clr <= 1e-3)) # check there are no 0s
## [1] 4
hseq_ctp_pheno.tmp.clr0 <- hseq_ctp_pheno.tmp.clr
hseq_ctp_pheno.tmp.clr0[which(hseq_ctp_pheno.tmp.clr0 <= 1e-3)] <- 1e-3
hseq_ctp_pheno.tmp.clr0.tr <- clr(hseq_ctp_pheno.tmp.clr0)
#hseq_ctp_pheno.clr.1e3 <- data.frame(keep_pheno_col, hseq_ctp_pheno.tmp.clr0.tr)
hseq_ctp_pheno.clr.1e3.tmp <- inner_join(data.frame(keep_pheno_col, hseq_ctp_pheno.tmp.clr0.tr), ctp_pheno[,c("IID", "projid")], by = "IID")
ctp_pheno_ihc.clr.1e3 <- inner_join(hseq_ctp_pheno.clr.1e3.tmp, ihc.clr.1e3, by = "projid")
ctp_pheno_ihc_hseq.clr.1e3 <- ctp_pheno_ihc.clr.1e3[,c("IID", all_of(c("hseq_Neuron", "hseq_Astro", "hseq_Endo", "hseq_Micro", "hseq_Oligo", "hseq_OPC")), "neuro_ihc", "astro_ihc", "endo_ihc", "micro_ihc", "oligo_ihc")]
ihc_rmoligo <- ihc %>%
mutate(neuro_ihc = neuro_ihc + oligo_ihc) %>%
dplyr::select(-c("oligo_ihc"))
# clr transform
ihc_rmoligo.clr.1e3.tmp0 <- as.matrix(ihc[,c("neuro_ihc", "astro_ihc", "endo_ihc", "micro_ihc")])
ihc_rmoligo.clr.1e3.tmp0[which(ihc_rmoligo.clr.1e3.tmp0 < 1e-3)] <- 1e-3
ihc_rmoligo.clr.1e3.tmp <- clr(ihc_rmoligo.clr.1e3.tmp0)
ihc_rmoligo.clr.1e3 <- data.frame(projid = ihc$projid, ihc_rmoligo.clr.1e3.tmp)
# Redo the hseq clr-transform with Neuron sum
hseq_ctp_pheno_adjoligo.tmp.clr <- as.matrix(hseq_ctp_pheno_adjoligo.tmp %>% dplyr::select(c("hseq_Neuron", "hseq_Astro", "hseq_Endo", "hseq_Micro", "hseq_OPC")))
length(which(hseq_ctp_pheno_adjoligo.tmp.clr <= 1e-3)) # check there are no 0s
## [1] 4
hseq_ctp_pheno_adjoligo.tmp.clr0 <- hseq_ctp_pheno_adjoligo.tmp.clr
hseq_ctp_pheno_adjoligo.tmp.clr0[which(hseq_ctp_pheno_adjoligo.tmp.clr0 <= 1e-3)] <- 1e-3
hseq_ctp_pheno_adjoligo.tmp.clr0.tr <- clr(hseq_ctp_pheno_adjoligo.tmp.clr0)
hseq_ctp_pheno_adjoligo.clr.1e3.tmp <- inner_join(data.frame(keep_pheno_col, hseq_ctp_pheno_adjoligo.tmp.clr0.tr), ctp_pheno[,c("IID", "projid")], by = "IID")
ctp_pheno_ihc_rmoligo.clr.1e3 <- inner_join(hseq_ctp_pheno_adjoligo.clr.1e3.tmp, ihc_rmoligo.clr.1e3, by = "projid")
ctp_pheno_ihc_rmoligo_hseq.clr.1e3 <- ctp_pheno_ihc_rmoligo.clr.1e3[,c("IID", all_of(c("hseq_Neuron", "hseq_Astro", "hseq_Endo", "hseq_Micro", "hseq_OPC")), "neuro_ihc", "astro_ihc", "endo_ihc", "micro_ihc")]
if (write_out == TRUE) {
write.table(ctp_pheno_ihc, paste(out_dir, "/", filen, "_ctp_pheno_ihc.txt", sep = ""), col.names = T, row.names = F, quote = F, sep = "\t")
write.table(ctp_pheno_ihc_hseq, paste(out_dir, "/", filen, "_ctp_pheno_hseq_ihc.txt", sep = ""), col.names = T, row.names = F, quote = F, sep = "\t")
write.table(ctp_pheno_ihc_hseq_rmoligo, paste(out_dir, "/", filen, "_ctp_pheno_hseq_ihc_rmoligo.txt", sep = ""), col.names = T, row.names = F, quote = F, sep = "\t")
write.table(ctp_pheno_ihc_hseq.clr.1e3, paste(out_dir, "/", filen, "_ctp_pheno_hseq_ihc_clr.1e3.txt", sep = ""), col.names = T, row.names = F, quote = F, sep = "\t")
write.table(ctp_pheno_ihc_rmoligo_hseq.clr.1e3, paste(out_dir, "/", filen, "_ctp_pheno_hseq_ihc_rmoligo_clr.1e3.txt", sep = ""), col.names = T, row.names = F, quote = F, sep = "\t")
}
print("Neurons")
## [1] "Neurons"
cor.test(ctp_pheno_ihc$neuro_ihc, ctp_pheno_ihc$hseq_Neuron)
##
## Pearson's product-moment correlation
##
## data: ctp_pheno_ihc$neuro_ihc and ctp_pheno_ihc$hseq_Neuron
## t = -0.39395, df = 47, p-value = 0.6954
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3331899 0.2274983
## sample estimates:
## cor
## -0.0573682
cor.test(ctp_pheno_ihc$neuro_ihc, ctp_pheno_ihc$mcc_Neuron)
##
## Pearson's product-moment correlation
##
## data: ctp_pheno_ihc$neuro_ihc and ctp_pheno_ihc$mcc_Neuron
## t = 0.84458, df = 47, p-value = 0.4026
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1645857 0.3900554
## sample estimates:
## cor
## 0.1222697
cor.test(ctp_pheno_ihc$neuro_ihc, ctp_pheno_ihc$h_Neuron)
##
## Pearson's product-moment correlation
##
## data: ctp_pheno_ihc$neuro_ihc and ctp_pheno_ihc$h_Neuron
## t = -0.16098, df = 47, p-value = 0.8728
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3026729 0.2594349
## sample estimates:
## cor
## -0.02347421
print("Astrocytes")
## [1] "Astrocytes"
cor.test(ctp_pheno_ihc$astro_ihc, ctp_pheno_ihc$hseq_Astro)
##
## Pearson's product-moment correlation
##
## data: ctp_pheno_ihc$astro_ihc and ctp_pheno_ihc$hseq_Astro
## t = 1.7801, df = 47, p-value = 0.08152
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.03214714 0.49736851
## sample estimates:
## cor
## 0.2513213
cor.test(ctp_pheno_ihc$astro_ihc, ctp_pheno_ihc$mcc_Astro)
##
## Pearson's product-moment correlation
##
## data: ctp_pheno_ihc$astro_ihc and ctp_pheno_ihc$mcc_Astro
## t = 1.9765, df = 47, p-value = 0.05398
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.004531529 0.517874043
## sample estimates:
## cor
## 0.277018
cor.test(ctp_pheno_ihc$astro_ihc, ctp_pheno_ihc$h_ASTRO)
##
## Pearson's product-moment correlation
##
## data: ctp_pheno_ihc$astro_ihc and ctp_pheno_ihc$h_ASTRO
## t = 0.6339, df = 47, p-value = 0.5292
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1941526 0.3638472
## sample estimates:
## cor
## 0.09207058
print("Endothelial")
## [1] "Endothelial"
cor.test(ctp_pheno_ihc$endo_ihc, ctp_pheno_ihc$hseq_Endo)
##
## Pearson's product-moment correlation
##
## data: ctp_pheno_ihc$endo_ihc and ctp_pheno_ihc$hseq_Endo
## t = 1.4661, df = 47, p-value = 0.1493
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.07657583 0.46308851
## sample estimates:
## cor
## 0.2091238
cor.test(ctp_pheno_ihc$endo_ihc, ctp_pheno_ihc$mcc_Endo)
##
## Pearson's product-moment correlation
##
## data: ctp_pheno_ihc$endo_ihc and ctp_pheno_ihc$mcc_Endo
## t = 0.97835, df = 47, p-value = 0.3329
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1457092 0.4063304
## sample estimates:
## cor
## 0.1412758
cor.test(ctp_pheno_ihc$endo_ihc, ctp_pheno_ihc$h_ENDO)
##
## Pearson's product-moment correlation
##
## data: ctp_pheno_ihc$endo_ihc and ctp_pheno_ihc$h_ENDO
## t = 1.1408, df = 47, p-value = 0.2598
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1227221 0.4256877
## sample estimates:
## cor
## 0.1641388
print("Microglia")
## [1] "Microglia"
cor.test(ctp_pheno_ihc$micro_ihc, ctp_pheno_ihc$hseq_Micro)
##
## Pearson's product-moment correlation
##
## data: ctp_pheno_ihc$micro_ihc and ctp_pheno_ihc$hseq_Micro
## t = 2.3986, df = 47, p-value = 0.02048
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.05407108 0.55948672
## sample estimates:
## cor
## 0.3302464
cor.test(ctp_pheno_ihc$micro_ihc, ctp_pheno_ihc$mcc_Micro)
##
## Pearson's product-moment correlation
##
## data: ctp_pheno_ihc$micro_ihc and ctp_pheno_ihc$mcc_Micro
## t = 1.0132, df = 47, p-value = 0.3161
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1407753 0.4105274
## sample estimates:
## cor
## 0.146209
cor.test(ctp_pheno_ihc$micro_ihc, ctp_pheno_ihc$h_MONO)
##
## Pearson's product-moment correlation
##
## data: ctp_pheno_ihc$micro_ihc and ctp_pheno_ihc$h_MONO
## t = 2.0768, df = 47, p-value = 0.04331
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.009502971 0.528070036
## sample estimates:
## cor
## 0.2899248
print("Oligodendrocytes/OPCs")
## [1] "Oligodendrocytes/OPCs"
cor.test(ctp_pheno_ihc$oligo_ihc, ctp_pheno_ihc$hseq_Oligo)
##
## Pearson's product-moment correlation
##
## data: ctp_pheno_ihc$oligo_ihc and ctp_pheno_ihc$hseq_Oligo
## t = 0.63397, df = 47, p-value = 0.5292
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1941425 0.3638563
## sample estimates:
## cor
## 0.09208094
cor.test(ctp_pheno_ihc$oligo_ihc, ctp_pheno_ihc$mcc_Oligo)
##
## Pearson's product-moment correlation
##
## data: ctp_pheno_ihc$oligo_ihc and ctp_pheno_ihc$mcc_Oligo
## t = 0.44019, df = 47, p-value = 0.6618
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.2211046 0.3391614
## sample estimates:
## cor
## 0.06407582
cor.test(ctp_pheno_ihc$oligo_ihc, ctp_pheno_ihc$h_OLIGO)
##
## Pearson's product-moment correlation
##
## data: ctp_pheno_ihc$oligo_ihc and ctp_pheno_ihc$h_OLIGO
## t = -0.47013, df = 47, p-value = 0.6404
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3430129 0.2169551
## sample estimates:
## cor
## -0.06841526
cor.test(ctp_pheno_ihc$oligo_ihc, ctp_pheno_ihc$hseq_OPC)
##
## Pearson's product-moment correlation
##
## data: ctp_pheno_ihc$oligo_ihc and ctp_pheno_ihc$hseq_OPC
## t = -0.40909, df = 47, p-value = 0.6843
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3351492 0.2254058
## sample estimates:
## cor
## -0.05956635
cor.test(ctp_pheno_ihc$oligo_ihc, ctp_pheno_ihc$mcc_OPC)
##
## Pearson's product-moment correlation
##
## data: ctp_pheno_ihc$oligo_ihc and ctp_pheno_ihc$mcc_OPC
## t = -0.7638, df = 47, p-value = 0.4488
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3800883 0.1759488
## sample estimates:
## cor
## -0.1107261
cor.test(ctp_pheno_ihc$oligo_ihc, ctp_pheno_ihc$h_OPC)
##
## Pearson's product-moment correlation
##
## data: ctp_pheno_ihc$oligo_ihc and ctp_pheno_ihc$h_OPC
## t = 0.52129, df = 47, p-value = 0.6046
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.2098519 0.3495627
## sample estimates:
## cor
## 0.07581873
adj <- "noadjoligo"
transf <- "noclr"
ctp_ihc.tmp <- ctp_pheno_ihc
neuron_cor <- cor.test(ctp_ihc.tmp$neuro_ihc, ctp_ihc.tmp$hseq_Neuron)
astro_cor <- cor.test(ctp_ihc.tmp$astro_ihc, ctp_ihc.tmp$hseq_Astro)
endo_cor <- cor.test(ctp_ihc.tmp$endo_ihc, ctp_ihc.tmp$hseq_Endo)
micro_cor <- cor.test(ctp_ihc.tmp$micro_ihc, ctp_ihc.tmp$hseq_Micro)
oligo_cor <- cor.test(ctp_ihc.tmp$oligo_ihc, ctp_ihc.tmp$hseq_Oligo)
neuro_ihc.long <- melt(ctp_ihc.tmp, measure.vars = c("hseq_Neuron", "mcc_Neuron", "h_Neuron", "h_NeuN_pos"), variable.name = "method", value.name = "CTP")
astro_ihc.long <- melt(ctp_ihc.tmp, measure.vars = c("hseq_Astro", "mcc_Astro", "h_ASTRO", "h_NeuN_neg"), variable.name = "method", value.name = "CTP")
endo_ihc.long <- melt(ctp_ihc.tmp, measure.vars = c("hseq_Endo", "mcc_Endo", "h_ENDO", "h_NeuN_neg"), variable.name = "method", value.name = "CTP")
micro_ihc.long <- melt(ctp_ihc.tmp, measure.vars = c("hseq_Micro", "mcc_Micro", "h_MONO", "h_NeuN_neg"), variable.name = "method", value.name = "CTP")
oligo_ihc.long <- melt(ctp_ihc.tmp, measure.vars = c("hseq_Oligo", "mcc_Oligo", "h_OLIGO", "h_NeuN_neg"), variable.name = "method", value.name = "CTP")
ggplot(neuro_ihc.long, aes(x = neuro_ihc, y = CTP)) + geom_point() + geom_smooth(method = "lm") + facet_wrap(~ method, scales = "free_y")
## `geom_smooth()` using formula 'y ~ x'
ggplot(astro_ihc.long, aes(x = astro_ihc, y = CTP)) + geom_point() + geom_smooth(method = "lm") + facet_wrap(~ method, scales = "free_y")
## `geom_smooth()` using formula 'y ~ x'
ggplot(endo_ihc.long, aes(x = endo_ihc, y = CTP)) + geom_point() + geom_smooth(method = "lm") + facet_wrap(~ method, scales = "free_y")
## `geom_smooth()` using formula 'y ~ x'
ggplot(micro_ihc.long, aes(x = micro_ihc, y = CTP)) + geom_point() + geom_smooth(method = "lm") + facet_wrap(~ method, scales = "free_y")
## `geom_smooth()` using formula 'y ~ x'
ggplot(oligo_ihc.long, aes(x = oligo_ihc, y = CTP)) + geom_point() + geom_smooth(method = "lm") + facet_wrap(~ method, scales = "free_y")
## `geom_smooth()` using formula 'y ~ x'
# Put hseq correlations into 1 plot
neuro_ihc.gg <- ggplot(neuro_ihc.long %>% filter(method == "hseq_Neuron"), aes(x = neuro_ihc, y = CTP)) + geom_point() + geom_smooth(method = "lm") + ylab("Neuron") + ggtitle(paste("r = ", signif(neuron_cor$estimate, 2), ", p = ", signif(neuron_cor$p.value, 2), sep = ""))
astro_ihc.gg <- ggplot(astro_ihc.long %>% filter(method == "hseq_Astro"), aes(x = astro_ihc, y = CTP)) + geom_point() + geom_smooth(method = "lm") + ylab("Astro") + ggtitle(paste("r = ", signif(astro_cor$estimate, 2), ", p = ", signif(astro_cor$p.value, 2), sep = ""))
endo_ihc.gg <- ggplot(endo_ihc.long %>% filter(method == "hseq_Endo"), aes(x = endo_ihc, y = CTP)) + geom_point() + geom_smooth(method = "lm") + ylab("Endo") + ggtitle(paste("r = ", signif(endo_cor$estimate, 2), ", p = ", signif(endo_cor$p.value, 2), sep = ""))
micro_ihc.gg <- ggplot(micro_ihc.long %>% filter(method == "hseq_Micro"), aes(x = micro_ihc, y = CTP)) + geom_point() + geom_smooth(method = "lm") + ylab("Micro") + ggtitle(paste("r = ", signif(micro_cor$estimate, 2), ", p = ", signif(micro_cor$p.value, 2), sep = ""))
oligo_ihc.gg <- ggplot(oligo_ihc.long %>% filter(method == "hseq_Oligo"), aes(x = oligo_ihc, y = CTP)) + geom_point() + geom_smooth(method = "lm") + ylab("Oligo") + ggtitle(paste("r = ", signif(oligo_cor$estimate, 2), ", p = ", signif(oligo_cor$p.value, 2), sep = ""))
ggarrange(neuro_ihc.gg, astro_ihc.gg, endo_ihc.gg, micro_ihc.gg, oligo_ihc.gg, ncol = 3, nrow = 2)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
ggsave(paste(out_dir, "/ctp_assoc_ihc_", adj, "_", transf, ".png", sep = ""), bg = "transparent", height = 7, width = 7)
adj <- "adjoligo"
transf <- "noclr"
ctp_ihc.tmp <- ctp_pheno_ihc_hseq_rmoligo
neuron_cor <- cor.test(ctp_ihc.tmp$neuro_ihc, ctp_ihc.tmp$hseq_Neuron)
astro_cor <- cor.test(ctp_ihc.tmp$astro_ihc, ctp_ihc.tmp$hseq_Astro)
endo_cor <- cor.test(ctp_ihc.tmp$endo_ihc, ctp_ihc.tmp$hseq_Endo)
micro_cor <- cor.test(ctp_ihc.tmp$micro_ihc, ctp_ihc.tmp$hseq_Micro)
neuro_ihc.long <- melt(ctp_ihc.tmp, measure.vars = c("hseq_Neuron"), variable.name = "method", value.name = "CTP")
astro_ihc.long <- melt(ctp_ihc.tmp, measure.vars = c("hseq_Astro"), variable.name = "method", value.name = "CTP")
endo_ihc.long <- melt(ctp_ihc.tmp, measure.vars = c("hseq_Endo"), variable.name = "method", value.name = "CTP")
micro_ihc.long <- melt(ctp_ihc.tmp, measure.vars = c("hseq_Micro"), variable.name = "method", value.name = "CTP")
# Put hseq correlations into 1 plot
neuro_ihc.gg <- ggplot(neuro_ihc.long %>% filter(method == "hseq_Neuron"), aes(x = neuro_ihc, y = CTP)) + geom_point() + geom_smooth(method = "lm") + ylab("Neuron") + ggtitle(paste("r = ", signif(neuron_cor$estimate, 2), ", p = ", signif(neuron_cor$p.value, 2), sep = ""))
astro_ihc.gg <- ggplot(astro_ihc.long %>% filter(method == "hseq_Astro"), aes(x = astro_ihc, y = CTP)) + geom_point() + geom_smooth(method = "lm") + ylab("Astro") + ggtitle(paste("r = ", signif(astro_cor$estimate, 2), ", p = ", signif(astro_cor$p.value, 2), sep = ""))
endo_ihc.gg <- ggplot(endo_ihc.long %>% filter(method == "hseq_Endo"), aes(x = endo_ihc, y = CTP)) + geom_point() + geom_smooth(method = "lm") + ylab("Endo") + ggtitle(paste("r = ", signif(endo_cor$estimate, 2), ", p = ", signif(endo_cor$p.value, 2), sep = ""))
micro_ihc.gg <- ggplot(micro_ihc.long %>% filter(method == "hseq_Micro"), aes(x = micro_ihc, y = CTP)) + geom_point() + geom_smooth(method = "lm") + ylab("Micro") + ggtitle(paste("r = ", signif(micro_cor$estimate, 2), ", p = ", signif(micro_cor$p.value, 2), sep = ""))
ggarrange(neuro_ihc.gg, astro_ihc.gg, endo_ihc.gg, micro_ihc.gg)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
ggsave(paste(out_dir, "/ctp_assoc_ihc_", adj, "_", transf, ".png", sep = ""), bg = "transparent", height = 7, width = 7)
adj <- "noadjoligo"
transf <- "clr"
ctp_ihc.tmp <- ctp_pheno_ihc_hseq.clr.1e3
neuron_cor <- cor.test(ctp_ihc.tmp$neuro_ihc, ctp_ihc.tmp$hseq_Neuron)
astro_cor <- cor.test(ctp_ihc.tmp$astro_ihc, ctp_ihc.tmp$hseq_Astro)
endo_cor <- cor.test(ctp_ihc.tmp$endo_ihc, ctp_ihc.tmp$hseq_Endo)
micro_cor <- cor.test(ctp_ihc.tmp$micro_ihc, ctp_ihc.tmp$hseq_Micro)
oligo_cor <- cor.test(ctp_ihc.tmp$oligo_ihc, ctp_ihc.tmp$hseq_Oligo)
neuro_ihc.long <- melt(ctp_ihc.tmp, measure.vars = c("hseq_Neuron"), variable.name = "method", value.name = "CTP")
astro_ihc.long <- melt(ctp_ihc.tmp, measure.vars = c("hseq_Astro"), variable.name = "method", value.name = "CTP")
endo_ihc.long <- melt(ctp_ihc.tmp, measure.vars = c("hseq_Endo"), variable.name = "method", value.name = "CTP")
micro_ihc.long <- melt(ctp_ihc.tmp, measure.vars = c("hseq_Micro"), variable.name = "method", value.name = "CTP")
oligo_ihc.long <- melt(ctp_ihc.tmp, measure.vars = c("hseq_Oligo"), variable.name = "method", value.name = "CTP")
# Put hseq correlations into 1 plot
neuro_ihc.gg <- ggplot(neuro_ihc.long %>% filter(method == "hseq_Neuron"), aes(x = neuro_ihc, y = CTP)) + geom_point() + geom_smooth(method = "lm") + ylab("Neuron") + ggtitle(paste("r = ", signif(neuron_cor$estimate, 2), ", p = ", signif(neuron_cor$p.value, 2), sep = ""))
astro_ihc.gg <- ggplot(astro_ihc.long %>% filter(method == "hseq_Astro"), aes(x = astro_ihc, y = CTP)) + geom_point() + geom_smooth(method = "lm") + ylab("Astro") + ggtitle(paste("r = ", signif(astro_cor$estimate, 2), ", p = ", signif(astro_cor$p.value, 2), sep = ""))
endo_ihc.gg <- ggplot(endo_ihc.long %>% filter(method == "hseq_Endo"), aes(x = endo_ihc, y = CTP)) + geom_point() + geom_smooth(method = "lm") + ylab("Endo") + ggtitle(paste("r = ", signif(endo_cor$estimate, 2), ", p = ", signif(endo_cor$p.value, 2), sep = ""))
micro_ihc.gg <- ggplot(micro_ihc.long %>% filter(method == "hseq_Micro"), aes(x = micro_ihc, y = CTP)) + geom_point() + geom_smooth(method = "lm") + ylab("Micro") + ggtitle(paste("r = ", signif(micro_cor$estimate, 2), ", p = ", signif(micro_cor$p.value, 2), sep = ""))
oligo_ihc.gg <- ggplot(oligo_ihc.long %>% filter(method == "hseq_Oligo"), aes(x = oligo_ihc, y = CTP)) + geom_point() + geom_smooth(method = "lm") + ylab("Oligo") + ggtitle(paste("r = ", signif(oligo_cor$estimate, 2), ", p = ", signif(oligo_cor$p.value, 2), sep = ""))
ggarrange(neuro_ihc.gg, astro_ihc.gg, endo_ihc.gg, micro_ihc.gg, oligo_ihc.gg, ncol = 3, nrow = 2)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
ggsave(paste(out_dir, "/ctp_assoc_ihc_", adj, "_", transf, ".png", sep = ""), bg = "transparent", height = 7, width = 7)
adj <- "adjoligo"
transf <- "clr"
ctp_ihc.tmp <- ctp_pheno_ihc_rmoligo_hseq.clr.1e3
neuron_cor <- cor.test(ctp_ihc.tmp$neuro_ihc, ctp_ihc.tmp$hseq_Neuron)
astro_cor <- cor.test(ctp_ihc.tmp$astro_ihc, ctp_ihc.tmp$hseq_Astro)
endo_cor <- cor.test(ctp_ihc.tmp$endo_ihc, ctp_ihc.tmp$hseq_Endo)
micro_cor <- cor.test(ctp_ihc.tmp$micro_ihc, ctp_ihc.tmp$hseq_Micro)
neuro_ihc.long <- melt(ctp_ihc.tmp, measure.vars = c("hseq_Neuron"), variable.name = "method", value.name = "CTP")
astro_ihc.long <- melt(ctp_ihc.tmp, measure.vars = c("hseq_Astro"), variable.name = "method", value.name = "CTP")
endo_ihc.long <- melt(ctp_ihc.tmp, measure.vars = c("hseq_Endo"), variable.name = "method", value.name = "CTP")
micro_ihc.long <- melt(ctp_ihc.tmp, measure.vars = c("hseq_Micro"), variable.name = "method", value.name = "CTP")
# Put hseq correlations into 1 plot
neuro_ihc.gg <- ggplot(neuro_ihc.long %>% filter(method == "hseq_Neuron"), aes(x = neuro_ihc, y = CTP)) + geom_point() + geom_smooth(method = "lm") + ylab("Neuron") + ggtitle(paste("r = ", signif(neuron_cor$estimate, 2), ", p = ", signif(neuron_cor$p.value, 2), sep = ""))
astro_ihc.gg <- ggplot(astro_ihc.long %>% filter(method == "hseq_Astro"), aes(x = astro_ihc, y = CTP)) + geom_point() + geom_smooth(method = "lm") + ylab("Astro") + ggtitle(paste("r = ", signif(astro_cor$estimate, 2), ", p = ", signif(astro_cor$p.value, 2), sep = ""))
endo_ihc.gg <- ggplot(endo_ihc.long %>% filter(method == "hseq_Endo"), aes(x = endo_ihc, y = CTP)) + geom_point() + geom_smooth(method = "lm") + ylab("Endo") + ggtitle(paste("r = ", signif(endo_cor$estimate, 2), ", p = ", signif(endo_cor$p.value, 2), sep = ""))
micro_ihc.gg <- ggplot(micro_ihc.long %>% filter(method == "hseq_Micro"), aes(x = micro_ihc, y = CTP)) + geom_point() + geom_smooth(method = "lm") + ylab("Micro") + ggtitle(paste("r = ", signif(micro_cor$estimate, 2), ", p = ", signif(micro_cor$p.value, 2), sep = ""))
ggarrange(neuro_ihc.gg, astro_ihc.gg, endo_ihc.gg, micro_ihc.gg)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
ggsave(paste(out_dir, "/ctp_assoc_ihc_", adj, "_", transf, ".png", sep = ""), bg = "transparent", height = 7, width = 7)
#ctp_pheno_ihc_mcc_rmoligo <- ctp_pheno_ihc_mcc %>% dplyr::select(-c("NonN_Oligo_MBP"))
# Long format for deconvolved proportions
ctp_pheno_ihc_hseq_rmoligo.long <- melt(ctp_pheno_ihc_hseq_rmoligo, measure.vars = c("hseq_Neuron", "hseq_Exc", "hseq_Inh", "hseq_Astro", "hseq_Endo", "hseq_Micro", "hseq_OPC"), variable.name = "celltype", value.name = "CTP_adjoligo")
# Long format for IHC proportions
ctp_pheno_ihc_hseq_rmoligo.long <- melt(ctp_pheno_ihc_hseq_rmoligo.long, measure.vars = c("neuro_ihc", "astro_ihc", "endo_ihc", "micro_ihc"), variable.name = "celltype_ihc", value.name = "CTP_ihc_adjoligo")
# Only take matching data
ctp_pheno_ihc_hseq_rmoligo.long <- ctp_pheno_ihc_hseq_rmoligo.long %>% filter((celltype == "hseq_Neuron" & celltype_ihc == "neuro_ihc") | (celltype == "hseq_Astro" & celltype_ihc == "astro_ihc") | (celltype == "hseq_Endo" & celltype_ihc == "endo_ihc") | (celltype == "hseq_Micro" & celltype_ihc == "micro_ihc"))
neuro_ihc_rmoligo <- ctp_pheno_ihc_hseq_rmoligo.long %>% filter(celltype_ihc == "neuro_ihc") %>% dplyr::select(c("CTP_adjoligo", "CTP_ihc_adjoligo"))
cor.test(neuro_ihc_rmoligo$CTP_ihc_adjoligo, neuro_ihc_rmoligo$CTP_adjoligo)
astro_ihc_rmoligo <- ctp_pheno_ihc_hseq_rmoligo.long %>% filter(celltype_ihc == "astro_ihc") %>% dplyr::select(c("CTP_adjoligo", "CTP_ihc_adjoligo"))
cor.test(astro_ihc_rmoligo$CTP_ihc_adjoligo, astro_ihc_rmoligo$CTP_adjoligo)
micro_ihc_rmoligo <- ctp_pheno_ihc_hseq_rmoligo.long %>% filter(celltype_ihc == "micro_ihc") %>% dplyr::select(c("CTP_adjoligo", "CTP_ihc_adjoligo"))
cor.test(micro_ihc_rmoligo$CTP_ihc_adjoligo, micro_ihc_rmoligo$CTP_adjoligo)
endo_ihc_rmoligo <- ctp_pheno_ihc_hseq_rmoligo.long %>% filter(celltype_ihc == "endo_ihc") %>% dplyr::select(c("CTP_adjoligo", "CTP_ihc_adjoligo"))
cor.test(endo_ihc_rmoligo$CTP_ihc_adjoligo, endo_ihc_rmoligo$CTP_adjoligo)
neuro_adjoligo.gg <- ggplot(ctp_pheno_ihc_hseq_rmoligo.long %>% filter(celltype_ihc == "neuro_ihc"), aes(x = CTP_ihc_adjoligo, y = CTP_adjoligo)) + geom_point() + geom_smooth(method = "lm") + ggtitle("neurons, adjusting for oligos")
astro_adjoligo.gg <- ggplot(ctp_pheno_ihc_hseq_rmoligo.long %>% filter(celltype_ihc == "astro_ihc"), aes(x = CTP_ihc_adjoligo, y = CTP_adjoligo)) + geom_point() + geom_smooth(method = "lm") + ggtitle("astrocytes, adjusting for oligos")
micro_adjoligo.gg <- ggplot(ctp_pheno_ihc_hseq_rmoligo.long %>% filter(celltype_ihc == "micro_ihc"), aes(x = CTP_ihc_adjoligo, y = CTP_adjoligo)) + geom_point() + geom_smooth(method = "lm") + ggtitle("microglia, adjusting for oligos")
endo_adjoligo.gg <- ggplot(ctp_pheno_ihc_hseq_rmoligo.long %>% filter(celltype_ihc == "endo_ihc"), aes(x = CTP_ihc_adjoligo, y = CTP_adjoligo)) + geom_point() + geom_smooth(method = "lm") + ggtitle("endothelial, adjusting for oligos")
ggarrange(neuro_adjoligo.gg, astro_adjoligo.gg, micro_adjoligo.gg, endo_adjoligo.gg, nrow = 2, ncol = 2)
ctp_pheno_ihc_hseq_rmoligo_pheno <- inner_join(ctp_pheno_ihc_hseq_rmoligo, pheno, by = "IID")
table(ctp_pheno_ihc_hseq_rmoligo_pheno$group)
##
## AZD CTL
## 18 31
keep_pheno_col <- ctp_pheno_ihc_hseq_rmoligo_pheno %>% dplyr::select(all_of(keep_cols))
ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr <- as.matrix(ctp_pheno_ihc_hseq_rmoligo_pheno %>% dplyr::select(c("neuro_ihc", "astro_ihc", "endo_ihc", "micro_ihc")))
# Make minimum CTP = 1e-3
length(which(ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr <= 1e-3)) # check there are no 0s
## [1] 0
ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr0 <- ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr
ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr0[which(ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr0 <= 1e-3)] <- 1e-3
# Perform clr-transform
ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr0.tr <- clr(ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr0)
ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3 <- data.frame(keep_pheno_col, ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr0.tr)
summary(lm(neuro_ihc ~ group, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3))
##
## Call:
## lm(formula = neuro_ihc ~ group, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.44023 -0.10196 0.03129 0.08372 0.34713
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.07393 0.03902 27.526 <2e-16 ***
## groupCTL 0.05423 0.04905 1.105 0.275
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1655 on 47 degrees of freedom
## Multiple R-squared: 0.02534, Adjusted R-squared: 0.004606
## F-statistic: 1.222 on 1 and 47 DF, p-value: 0.2746
summary(lm(astro_ihc ~ group, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3))
##
## Call:
## lm(formula = astro_ihc ~ group, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.32685 -0.14919 0.01582 0.15250 0.32810
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.11217 0.04211 -2.664 0.0106 *
## groupCTL -0.11212 0.05295 -2.118 0.0395 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1787 on 47 degrees of freedom
## Multiple R-squared: 0.08709, Adjusted R-squared: 0.06767
## F-statistic: 4.484 on 1 and 47 DF, p-value: 0.03953
summary(lm(endo_ihc ~ group, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3))
##
## Call:
## lm(formula = endo_ihc ~ group, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.23531 -0.13337 -0.00206 0.09637 0.42833
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.17217 0.03708 -4.643 2.78e-05 ***
## groupCTL 0.05204 0.04662 1.116 0.27
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1573 on 47 degrees of freedom
## Multiple R-squared: 0.02582, Adjusted R-squared: 0.005095
## F-statistic: 1.246 on 1 and 47 DF, p-value: 0.27
summary(lm(micro_ihc ~ group, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3))
##
## Call:
## lm(formula = micro_ihc ~ group, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.48884 -0.15464 0.00626 0.13685 0.54099
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.789579 0.052818 -14.949 <2e-16 ***
## groupCTL 0.005854 0.066404 0.088 0.93
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2241 on 47 degrees of freedom
## Multiple R-squared: 0.0001653, Adjusted R-squared: -0.02111
## F-statistic: 0.007771 on 1 and 47 DF, p-value: 0.9301
summary(lm(neuro_ihc ~ group + age + batch + race + pmi, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3))
##
## Call:
## lm(formula = neuro_ihc ~ group + age + batch + race + pmi, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.41434 -0.10129 0.03912 0.10099 0.31929
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.620323 0.490546 1.265 0.213
## groupCTL 0.083826 0.053849 1.557 0.127
## age 0.004322 0.005358 0.807 0.424
## batch1 -0.012311 0.050227 -0.245 0.808
## race 0.101193 0.103188 0.981 0.332
## pmi -0.005057 0.005537 -0.913 0.366
##
## Residual standard error: 0.1684 on 43 degrees of freedom
## Multiple R-squared: 0.07748, Adjusted R-squared: -0.02979
## F-statistic: 0.7223 on 5 and 43 DF, p-value: 0.6104
summary(lm(astro_ihc ~ group + age + batch + race + pmi, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3))
##
## Call:
## lm(formula = astro_ihc ~ group + age + batch + race + pmi, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.30858 -0.14951 0.01124 0.14116 0.34214
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8.137e-02 5.399e-01 -0.151 0.8809
## groupCTL -1.144e-01 5.927e-02 -1.930 0.0603 .
## age 1.784e-05 5.897e-03 0.003 0.9976
## batch1 -4.490e-02 5.529e-02 -0.812 0.4212
## race -5.564e-03 1.136e-01 -0.049 0.9612
## pmi 2.467e-04 6.095e-03 0.040 0.9679
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1853 on 43 degrees of freedom
## Multiple R-squared: 0.1016, Adjusted R-squared: -0.002912
## F-statistic: 0.9721 on 5 and 43 DF, p-value: 0.4456
summary(lm(endo_ihc ~ group + age + batch + race + pmi, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3))
##
## Call:
## lm(formula = endo_ihc ~ group + age + batch + race + pmi, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.26012 -0.11277 0.00086 0.08934 0.43386
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.130535 0.473790 0.276 0.784
## groupCTL 0.051419 0.052009 0.989 0.328
## age -0.003147 0.005175 -0.608 0.546
## batch1 -0.004344 0.048512 -0.090 0.929
## race -0.004672 0.099663 -0.047 0.963
## pmi -0.004291 0.005348 -0.802 0.427
##
## Residual standard error: 0.1626 on 43 degrees of freedom
## Multiple R-squared: 0.04796, Adjusted R-squared: -0.06275
## F-statistic: 0.4332 on 5 and 43 DF, p-value: 0.8229
summary(lm(micro_ihc ~ group + age + batch + race + pmi, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3))
##
## Call:
## lm(formula = micro_ihc ~ group + age + batch + race + pmi, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.42155 -0.18108 0.00904 0.12077 0.49763
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.669486 0.661677 -1.012 0.317
## groupCTL -0.020876 0.072634 -0.287 0.775
## age -0.001193 0.007227 -0.165 0.870
## batch1 0.061551 0.067749 0.909 0.369
## race -0.090957 0.139186 -0.653 0.517
## pmi 0.009101 0.007469 1.219 0.230
##
## Residual standard error: 0.2271 on 43 degrees of freedom
## Multiple R-squared: 0.0605, Adjusted R-squared: -0.04875
## F-statistic: 0.5538 on 5 and 43 DF, p-value: 0.7346
ctp_pheno_ihc_hseq_rmoligo_pheno <- inner_join(ctp_pheno_ihc_hseq_rmoligo, pheno, by = "IID")
table(ctp_pheno_ihc_hseq_rmoligo_pheno$group)
##
## AZD CTL
## 18 31
keep_pheno_col <- ctp_pheno_ihc_hseq_rmoligo_pheno %>% dplyr::select(all_of(keep_cols))
ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr <- as.matrix(ctp_pheno_ihc_hseq_rmoligo_pheno %>% dplyr::select(c("hseq_Exc", "hseq_Inh", "hseq_Astro", "hseq_Endo", "hseq_Micro", "hseq_OPC")))
# Make minimum CTP = 1e-3
length(which(ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr <= 1e-3)) # check there are no 0s
## [1] 0
ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr0 <- ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr
ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr0[which(ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr0 <= 1e-3)] <- 1e-3
# Perform clr-transform
ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr0.tr <- clr(ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr0)
ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3 <- data.frame(keep_pheno_col, ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr0.tr)
summary(lm(hseq_Exc ~ group, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3))
##
## Call:
## lm(formula = hseq_Exc ~ group, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.34293 -0.08552 0.01162 0.07871 0.56493
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.30559 0.04196 31.115 <2e-16 ***
## groupCTL 0.03805 0.05275 0.721 0.474
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.178 on 47 degrees of freedom
## Multiple R-squared: 0.01095, Adjusted R-squared: -0.0101
## F-statistic: 0.5202 on 1 and 47 DF, p-value: 0.4743
summary(lm(hseq_Inh ~ group, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3))
##
## Call:
## lm(formula = hseq_Inh ~ group, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.133640 -0.054484 -0.008477 0.053479 0.290462
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.01618 0.02059 49.344 <2e-16 ***
## groupCTL 0.01731 0.02589 0.668 0.507
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08737 on 47 degrees of freedom
## Multiple R-squared: 0.009416, Adjusted R-squared: -0.01166
## F-statistic: 0.4467 on 1 and 47 DF, p-value: 0.5072
summary(lm(hseq_Astro ~ group, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3))
##
## Call:
## lm(formula = hseq_Astro ~ group, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.25131 -0.07718 -0.03181 0.09109 0.47219
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.07528 0.03386 31.753 <2e-16 ***
## groupCTL -0.02187 0.04258 -0.514 0.61
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1437 on 47 degrees of freedom
## Multiple R-squared: 0.005581, Adjusted R-squared: -0.01558
## F-statistic: 0.2638 on 1 and 47 DF, p-value: 0.6099
summary(lm(hseq_Endo ~ group, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3))
##
## Call:
## lm(formula = hseq_Endo ~ group, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.168696 -0.059976 -0.002166 0.053907 0.270680
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.64656 0.01897 -34.082 <2e-16 ***
## groupCTL -0.02342 0.02385 -0.982 0.331
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08049 on 47 degrees of freedom
## Multiple R-squared: 0.02011, Adjusted R-squared: -0.0007413
## F-statistic: 0.9644 on 1 and 47 DF, p-value: 0.3311
summary(lm(hseq_Micro ~ group, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3))
##
## Call:
## lm(formula = hseq_Micro ~ group, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.39353 -0.08995 0.00467 0.10263 0.24708
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.873426 0.032502 -26.873 <2e-16 ***
## groupCTL -0.004447 0.040862 -0.109 0.914
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1379 on 47 degrees of freedom
## Multiple R-squared: 0.000252, Adjusted R-squared: -0.02102
## F-statistic: 0.01185 on 1 and 47 DF, p-value: 0.9138
summary(lm(hseq_Exc ~ group + age + batch + race + pmi, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3))
##
## Call:
## lm(formula = hseq_Exc ~ group + age + batch + race + pmi, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.39429 -0.06578 0.01748 0.08145 0.52352
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.997912 0.496730 2.009 0.0508 .
## groupCTL 0.046405 0.054527 0.851 0.3995
## age 0.003187 0.005425 0.588 0.5599
## batch1 -0.129686 0.050860 -2.550 0.0144 *
## race 0.091031 0.104489 0.871 0.3885
## pmi 0.001836 0.005607 0.327 0.7449
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1705 on 43 degrees of freedom
## Multiple R-squared: 0.1701, Adjusted R-squared: 0.07362
## F-statistic: 1.763 on 5 and 43 DF, p-value: 0.141
summary(lm(hseq_Inh ~ group + age + batch + race + pmi, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3))
##
## Call:
## lm(formula = hseq_Inh ~ group + age + batch + race + pmi, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.141550 -0.069167 -0.000423 0.044592 0.273307
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.1469562 0.2594565 4.421 6.58e-05 ***
## groupCTL 0.0082869 0.0284813 0.291 0.772
## age -0.0014492 0.0028337 -0.511 0.612
## batch1 -0.0270007 0.0265658 -1.016 0.315
## race 0.0004041 0.0545774 0.007 0.994
## pmi 0.0021231 0.0029288 0.725 0.472
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08905 on 43 degrees of freedom
## Multiple R-squared: 0.05858, Adjusted R-squared: -0.05088
## F-statistic: 0.5352 on 5 and 43 DF, p-value: 0.7484
summary(lm(hseq_Astro ~ group + age + batch + race + pmi, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3))
##
## Call:
## lm(formula = hseq_Astro ~ group + age + batch + race + pmi, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.34161 -0.05518 -0.00325 0.06722 0.38626
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.002190 0.382236 2.622 0.01204 *
## groupCTL -0.009987 0.041959 -0.238 0.81300
## age 0.001555 0.004175 0.372 0.71136
## batch1 -0.135362 0.039137 -3.459 0.00124 **
## race 0.038974 0.080404 0.485 0.63033
## pmi -0.004336 0.004315 -1.005 0.32053
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1312 on 43 degrees of freedom
## Multiple R-squared: 0.2415, Adjusted R-squared: 0.1533
## F-statistic: 2.738 on 5 and 43 DF, p-value: 0.03111
summary(lm(hseq_Endo ~ group + age + batch + race + pmi, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3))
##
## Call:
## lm(formula = hseq_Endo ~ group + age + batch + race + pmi, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.164800 -0.050172 0.008049 0.036296 0.205258
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.5023769 0.2165868 -2.320 0.0252 *
## groupCTL -0.0182668 0.0237753 -0.768 0.4465
## age -0.0007989 0.0023655 -0.338 0.7372
## batch1 -0.0639759 0.0221764 -2.885 0.0061 **
## race -0.0089890 0.0455597 -0.197 0.8445
## pmi -0.0050250 0.0024449 -2.055 0.0460 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07434 on 43 degrees of freedom
## Multiple R-squared: 0.2353, Adjusted R-squared: 0.1464
## F-statistic: 2.646 on 5 and 43 DF, p-value: 0.03586
summary(lm(hseq_Micro ~ group + age + batch + race + pmi, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3))
##
## Call:
## lm(formula = hseq_Micro ~ group + age + batch + race + pmi, data = ctp_pheno_ihc_hseq_rmoligo_pheno.tmp.clr.1e3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.34702 -0.08913 -0.00494 0.10977 0.25531
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.9628320 0.4085811 -2.357 0.0231 *
## groupCTL -0.0158817 0.0448511 -0.354 0.7250
## age 0.0004300 0.0044624 0.096 0.9237
## batch1 0.0268570 0.0418347 0.642 0.5243
## race 0.0007982 0.0859462 0.009 0.9926
## pmi 0.0067727 0.0046121 1.468 0.1493
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1402 on 43 degrees of freedom
## Multiple R-squared: 0.05404, Adjusted R-squared: -0.05596
## F-statistic: 0.4913 on 5 and 43 DF, p-value: 0.7809